code: def excitations(electrons, orbitals, delta_sz=0):
      r"""Generate single and double excitations from a Hartree-Fock reference state.
  
      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_{r \in \mathrm{occ} \\ p \in \mathrm{unocc}}
          \hat{c}_p^\dagger \hat{c}_r \\
          && \hat{T}_2 = \sum_{r>s \in \mathrm{occ} \\ p>q \in
          \mathrm{unocc}} \hat{c}_p^\dagger \hat{c}_q^\dagger \hat{c}_r \hat{c}_s.
  
  
      In the equations above the indices :math:`r, s` and :math:`p, q` 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.
  
      |
  
      .. figure:: ../../_static/qchem/sd_excitations.png
          :align: center
          :width: 80%
  
      |
  
      Args:
          electrons (int): Number of electrons. If an active space is defined, this
              is the number of active electrons.
          orbitals (int): Number of *spin* orbitals. If an active space is defined,
              this is the number of active spin-orbitals.
          delta_sz (int): Specifies the selection rules ``sz[p] - sz[r] = delta_sz`` and
              ``sz[p] + sz[p] - sz[r] - sz[s] = delta_sz`` for the spin-projection ``sz`` of
              the orbitals involved in the single and double excitations, respectively.
              ``delta_sz`` can take the values :math:`0`, :math:`\pm 1` and :math:`\pm 2`.
  
      Returns:
          tuple(list, list): lists with the indices of the spin orbitals involved in the
          single and double excitations
  
      **Example**
  
      >>> electrons = 2
      >>> orbitals = 4
      >>> singles, doubles = excitations(electrons, orbitals)
      >>> print(singles)
      [[0, 2], [1, 3]]
      >>> print(doubles)
      [[0, 1, 2, 3]]
      """
  
      if not electrons > 0:
          raise ValueError(
              f"The number of active electrons has to be greater than 0 \n"
              f"Got n_electrons = {electrons}"
          )
  
      if orbitals <= electrons:
          raise ValueError(
              f"The number of active spin-orbitals ({orbitals}) "
              f"has to be greater than the number of active electrons ({electrons})."
          )
  
      if delta_sz not in (0, 1, -1, 2, -2):
          raise ValueError(
              f"Expected values for 'delta_sz' are 0, +/- 1 and +/- 2 but got ({delta_sz})."
          )
  
      # define the spin projection 'sz' of the single-particle states
      sz = np.array(
          [0.5 if (i % 2 == 0) else -0.5 for i in range(orbitals)]
      )
  
      singles = [
          [r, p]
          for r in range(electrons)
          for p in range(electrons, orbitals)
          if sz[p] - sz[r] == delta_sz
      ]
  
      doubles = [
          [s, r, q, p]
          for s in range(electrons - 1)
          for r in range(s + 1, electrons)
          for q in range(electrons, orbitals - 1)
          for p in range(q + 1, orbitals)
          if (sz[p] + sz[q] - sz[r] - sz[s]) == delta_sz
      ]
  
      return singles, doubles
  
file_path: pennylane/pennylane/qchem/structure.py
imports: import os
  import re
  from shutil import copyfile
  import numpy
label_desc: Use image node_id 3 for calling a global function with example usage: excitations(electrons, orbitals, delta_sz) and returns: singles,  doubles