Repository: PennyLaneAI/pennylane
Pull request #2227 (2022-02-24 10:30:09)

Make `grad_recipe` the first differentiation criterion in parameter shifts
This PR belongs to the series of #2180,  #2200, and #2207.

Here we make `grad_recipe` the first differentiation criterion to be checked in `param_shift`.
This allows users to easily define a hard-coded differentiation rule without having to pro-actively remove `parameter_frequencies` or the `generator`.
As noted in the new docstring for `qml.operation.Operation`, `grad_recipe` is not compatible with dynamic shift values and other context-based modifications of differentiation rules, so it should be used in the standard operations as little as possible.

As a side modification, the default two-term shift rule is removed from `Operation.get_parameter_shift`. Instead, an error is raised if no `grad_recipe` is available. `param_shift` does not make use of `get_parameter_shift` any longer, moving further towards its deprecation.

Related issues:
#2220 

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2227](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (e7c6cb1) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/1c0a57c68e8e7d237bbad29109053dd70e9d3716?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (1c0a57c) will **decrease** coverage by `2.75%`.
> The diff coverage is `95.71%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2227      +/-   ##
==========================================
- Coverage   99.26%   96.51%   -2.76%     
==========================================
  Files         231      231              
  Lines       18330    18310      -20     
==========================================
- Hits        18196    17671     -525     
- Misses        134      639     +505     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/ops/qubit/parametric\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9wYXJhbWV0cmljX29wcy5weQ==) | `100.00% <ø> (ø)` | |
| [...ylane/templates/subroutines/commuting\_evolution.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9jb21tdXRpbmdfZXZvbHV0aW9uLnB5) | `100.00% <ø> (ø)` | |
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `95.32% <90.90%> (-0.63%)` | :arrow_down: |
| [pennylane/tape/qubit\_param\_shift.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvcXViaXRfcGFyYW1fc2hpZnQucHk=) | `98.63% <91.66%> (-1.37%)` | :arrow_down: |
| [pennylane/gradients/parameter\_shift.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9wYXJhbWV0ZXJfc2hpZnQucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/cv.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9jdi5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/qubit/qchem\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9xY2hlbV9vcHMucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/tape/cv\_param\_shift.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvY3ZfcGFyYW1fc2hpZnQucHk=) | `100.00% <100.00%> (ø)` | |
| [...mplates/subroutines/fermionic\_double\_excitation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9mZXJtaW9uaWNfZG91YmxlX2V4Y2l0YXRpb24ucHk=) | `100.00% <100.00%> (ø)` | |
| [...mplates/subroutines/fermionic\_single\_excitation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9mZXJtaW9uaWNfc2luZ2xlX2V4Y2l0YXRpb24ucHk=) | `100.00% <100.00%> (ø)` | |
| ... and [24 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [1c0a57c...e7c6cb1](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2227?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 93a7c4c56..8b13d8b8f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -140,6 +140,7 @@
   rule is now automatically used by `qml.gradients.param_shift`.
   [(#2180)](https://github.com/PennyLaneAI/pennylane/pull/2180)
   [(#2182)](https://github.com/PennyLaneAI/pennylane/pull/2182)
+  [(#2227)](https://github.com/PennyLaneAI/pennylane/pull/2227)
 
   The frequencies can be used for circuit analysis, optimization
   via the `RotosolveOptimizer` and differentiation with the
@@ -150,14 +151,9 @@
   By default, the frequencies will be obtained from the
   `generator` property (if it is defined).
 
-  When using `qml.gradients.param_shift`, the parameter frequencies
-  are used to obtain the shift rule for the operation.
-
-  For operations that are registered to have an analytic gradient
-  method but that do not provide parameter frequencies, the
-  `grad_recipe` of the operation will be used for differentiation
-  instead. If there is no `grad_recipe`, the standard two-term shift
-  rule will be used.
+  When using `qml.gradients.param_shift`, either a custom `grad_recipe`
+  or the parameter frequencies are used to obtain the shift rule
+  for the operation, in that order of preference.
 
   See [Vidal and Theis (2018)](https://arxiv.org/abs/1812.06323)
   and [Wierichs et al. (2021)](https://arxiv.org/abs/2107.12390)
@@ -220,6 +216,15 @@
 
 <h3>Improvements</h3>
 
+* No two-term parameter-shift rule is assumed anymore by default.
+  [(#2227)](https://github.com/PennyLaneAI/pennylane/pull/2227)
+
+  Previously, operations marked for analytic differentiation that
+  do not provide a `generator`, `parameter_frequencies` or a
+  custom `grad_recipe` were assumed to satisfy the two-term shift
+  rule. This now has to be made explicit for custom operations
+  by adding any of the above attributes.
+
 * The `qml.draw_mpl` transform supports a `expansion_strategy` keyword argument.
   [(#2271)](https://github.com/PennyLaneAI/pennylane/pull/2271/)
 
@@ -253,6 +258,14 @@
 
 <h3>Deprecations</h3>
 
+* The `qml.operation.Operation.get_parameter_shift` method has been deprecated
+  and will be removed in a future release.
+  [#2227](https://github.com/PennyLaneAI/pennylane/pull/2227)
+
+  Instead, the functionalities for general parameter-shift rules in the
+  `qml.gradients` module should be used, together with the operation attributes
+  `parameter_frequencies` or `grad_recipe`.
+
 * The `qml.finite_diff()` function has been deprecated and will be removed
   in an upcoming release. Instead,
   `qml.gradients.finite_diff()` can be used to compute purely quantum gradients
diff --git a/pennylane/gradients/parameter_shift.py b/pennylane/gradients/parameter_shift.py
index 4cd1401dc..4c34bfbbb 100644
--- a/pennylane/gradients/parameter_shift.py
+++ b/pennylane/gradients/parameter_shift.py
@@ -73,42 +73,40 @@ def _get_operation_recipe(tape, t_idx, shifts):
 
     This function performs multiple attempts to obtain the recipe:
 
-    - If ``parameter_frequencies`` yield a result, the frequencies are
+    - If the operation has a custom :attr:`~.grad_recipe` defined, it is used.
+
+    - If :attr:`.parameter_frequencies` yields a result, the frequencies are
       used to construct the general parameter-shift rule via
-      ``qml.gradients.generate_shift_rule``
-
-    - If ``parameter_frequencies`` raises an error (because the operation
-      has no custom ``parameter_frequencies`` and no ``generator`` is defined),
-      the :meth:`.Operator.get_parameter_shift` method is used.
-      If in turn no custom version of ``get_parameter_shift`` is defined
-      and the operation does not have a ``grad_recipe``, the two-term
-      parameter-shift rule is assumed.
-
-    That is, a default to the two-term rule only is returned if no custom
-    ``parameter_frequencies``, ``generator``, ``get_parameter_shift``,
-    and ``grad_recipe`` are defined.
+      :func:`.generate_shift_rule`.
+      Note that by default, the generator is used to compute the parameter frequencies
+      if they are not provided by a custom implementation.
+
+    That is, the order of precedence is :meth:`~.grad_recipe`, custom
+    :attr:`~.parameter_frequencies`, and finally :meth:`.generator` via the default
+    implementation of the frequencies.
     """
     op, p_idx = tape.get_operation(t_idx)
+
+    # Try to use the stored grad_recipe of the operation
+    recipe = op.grad_recipe[p_idx]
+    if recipe is not None:
+        return process_shifts(np.array(recipe).T, check_duplicates=False)
+
+    # Try to obtain frequencies, either via custom implementation or from generator eigvals
     try:
-        # Obtain frequencies, either via custom implementation or from generator eigvals
         frequencies = op.parameter_frequencies[p_idx]
-        # Create shift rule from frequencies with given shifts
-        coeffs, shifts = qml.gradients.generate_shift_rule(frequencies, shifts=shifts, order=1)
-        # The shift rules do not include a rescaling of the parameter, only shifts.
-        mults = np.ones_like(coeffs)
-        return coeffs, mults, shifts
-
-    except qml.operation.OperatorPropertyUndefined:
-        # if no frequencies could be obtained, use the operation's shift rule if present.
-        # Note that the keyword argument ``shift`` only is used for the default two-term rule
-        if shifts is None:
-            use_shift = None
-        else:
-            use_shift = shifts[0]
+    except qml.operation.ParameterFrequenciesUndefinedError as e:
+        raise qml.operation.OperatorPropertyUndefined(
+            f"The operation {op.name} does not have a grad_recipe, parameter_frequencies or "
+            "a generator defined. No parameter shift rule can be applied."
+        ) from e
 
-        recipe = np.array(op.get_parameter_shift(p_idx, shift=use_shift)).T
+    # Create shift rule from frequencies with given shifts
+    coeffs, shifts = qml.gradients.generate_shift_rule(frequencies, shifts=shifts, order=1)
+    # The generated shift rules do not include a rescaling of the parameter, only shifts.
+    mults = np.ones_like(coeffs)
 
-        return process_shifts(recipe, check_duplicates=False)
+    return coeffs, mults, shifts
 
 
 def _gradient_analysis(tape, use_graph=True):
@@ -475,13 +473,21 @@ def param_shift(
         U(\mathbf{p})^\dagger \hat{O} U(\mathbf{p}) \vert 0\rangle.
 
 
-    The gradient of this expectation value can be calculated using :math:`2N` expectation
-    values using the parameter-shift rule:
+    The gradient of this expectation value can be calculated via the parameter-shift rule:
 
     .. math::
 
-        \frac{\partial f}{\partial \mathbf{p}} = \frac{1}{2\sin s} \left[ f(\mathbf{p} + s) -
-        f(\mathbf{p} -s) \right].
+        \frac{\partial f}{\partial \mathbf{p}} = \sum_{\mu=1}^{2R}
+        f\left(\mathbf{p}+\frac{2\mu-1}{2R}\pi\right)
+        \frac{(-1)^{\mu-1}}{4R\sin^2\left(\frac{2\mu-1}{4R}\pi\right)}
+
+    Here, :math:`R` is the number of frequencies with which the parameter :math:`\mathbf{p}`
+    enters the function :math:`f` via the operation :math:`U`, and we assumed that these
+    frequencies are equidistant.
+    For more general shift rules, both regarding the shifts and the frequencies, and
+    for more technical details, see
+    `Vidal and Theis (2018) <https://arxiv.org/abs/1812.06323>`_ and
+    `Wierichs et al. (2021) <https://arxiv.org/abs/2107.12390>`_.
 
     **Gradients of variances**
 
@@ -502,17 +508,16 @@ def param_shift(
         \mathbf{p}} \langle \hat{O}^2 \rangle (\mathbf{p})
         - 2 f(\mathbf{p}) \frac{\partial f}{\partial \mathbf{p}}.
 
-    This results in :math:`4N + 1` evaluations.
+    The derivatives in the expression on the right hand side can be computed via
+    the shift rule as above, allowing for the computation of the variance derivative.
 
-    In the case where :math:`O` is involutory (:math:`\hat{O}^2 = I`), the first term in the above
-    expression vanishes, and we are simply left with
+    In the case where :math:`O` is involutory (:math:`\hat{O}^2 = I`), the first
+    term in the above expression vanishes, and we are simply left with
 
     .. math::
 
       \frac{\partial g}{\partial \mathbf{p}} = - 2 f(\mathbf{p})
-      \frac{\partial f}{\partial \mathbf{p}},
-
-    allowing us to compute the gradient using :math:`2N + 1` evaluations.
+      \frac{\partial f}{\partial \mathbf{p}}.
 
     **Example**
 
@@ -531,6 +536,24 @@ def param_shift(
     tensor([[-0.38751725, -0.18884792, -0.38355708],
             [ 0.69916868,  0.34072432,  0.69202365]], requires_grad=True)
 
+    .. note::
+
+        ``param_shift`` performs multiple attempts to obtain the gradient recipes for
+        each operation:
+
+        - If an operation has a custom :attr:`~.operation.Operation.grad_recipe` defined,
+          it is used.
+
+        - If :attr:`~.operation.Operation.parameter_frequencies` yields a result, the frequencies
+          are used to construct the general parameter-shift rule via
+          :func:`.generate_shift_rule`.
+          Note that by default, the generator is used to compute the parameter frequencies
+          if they are not provided via a custom implementation.
+
+        That is, the order of precedence is :attr:`~.operation.Operation.grad_recipe`, custom
+        :attr:`~.operation.Operation.parameter_frequencies`, and finally
+        :meth:`~.operation.Operation.generator` via the default implementation of the frequencies.
+
     .. UsageDetails::
 
         This gradient transform can be applied directly to :class:`QNode <pennylane.QNode>` objects:
diff --git a/pennylane/operation.py b/pennylane/operation.py
index f86706324..487f0f88c 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -246,6 +246,11 @@ class GeneratorUndefinedError(OperatorPropertyUndefined):
     does not have a generator"""
 
 
+class ParameterFrequenciesUndefinedError(OperatorPropertyUndefined):
+    """Exception used to indicate that an operator
+    does not have parameter_frequencies"""
+
+
 # =============================================================================
 # Wire types
 # =============================================================================
@@ -1128,20 +1133,24 @@ class Operator(abc.ABC):
 class Operation(Operator):
     r"""Base class representing quantum gates or channels applied to quantum states.
 
-    Operations define some additional properties, such as differentiation
+    Operations define some additional properties, that are used for external
+    transformations such as gradient transforms.
 
-    The following two class attributes are optional, but in most cases
-    should be clearly defined to avoid unexpected behavior during
+    The following three class attributes are optional, but in most cases
+    at least one should be clearly defined to avoid unexpected behavior during
     differentiation.
 
-    * :attr:`~.Operation.grad_method`
     * :attr:`~.Operation.grad_recipe`
-
-    Finally, there are some additional optional class attributes
-    that may be set, and used by certain quantum optimizers:
-
+    * :attr:`~.Operation.parameter_frequencies`
     * :attr:`~.Operation.generator`
 
+    Note that ``grad_recipe`` takes precedence when computing parameter-shift
+    derivatives. Finally, these optional class attributes are used by certain
+    transforms, quantum optimizers, and gradient methods.
+    For details on how they are used during differentiation and other transforms,
+    please see the documentation for :class:`~.gradients.param_shift`,
+    :class:`~.metric_tensor`, :func:`~.reconstruct`.
+
     Args:
         params (tuple[tensor_like]): trainable parameters
         wires (Iterable[Any] or Any): Wire label(s) that the operator acts on.
@@ -1164,10 +1173,12 @@ class Operation(Operator):
         """
         if self.num_params == 0:
             return None
+        if self.grad_recipe != [None] * self.num_params:
+            return "A"
         try:
             self.parameter_frequencies  # pylint:disable=pointless-statement
             return "A"
-        except OperatorPropertyUndefined:
+        except ParameterFrequenciesUndefinedError:
             return "F"
 
     grad_recipe = None
@@ -1220,15 +1231,11 @@ class Operation(Operator):
         """
         raise NotImplementedError
 
-    def get_parameter_shift(self, idx, shift=None):
+    def get_parameter_shift(self, idx):
         r"""Multiplier and shift for the given parameter, based on its gradient recipe.
 
         Args:
-            idx (int): parameter index
-            shift (float or None): The shift value to use for the two-term
-                parameter-shift rule. This is only used if the operation
-                does not have :attr:`Operator.grad_recipe` defined. If ``None``,
-                a shift value of :math:`\pi/2` is used.
+            idx (int): parameter index within the operation
 
         Returns:
             list[[float, float, float]]: list of multiplier, coefficient, shift for each term in the gradient recipe
@@ -1236,25 +1243,22 @@ class Operation(Operator):
         Note that the default value for ``shift`` is None, which is replaced by the
         default shift :math:`\pi/2`.
         """
+        warnings.warn(
+            "The method get_parameter_shift is deprecated. Use the methods of "
+            "the gradients module for general parameter-shift rules instead.",
+            UserWarning,
+        )
         # get the gradient recipe for this parameter
         recipe = self.grad_recipe[idx]
         if recipe is not None:
             return recipe
 
-        # Default values
-        if shift is None:
-            shift = np.pi / 2
-        multiplier = 0.5 / np.sin(shift)
-        a = 1
-
-        # We set the following default recipe:
-        # ∂f(x) = c*f(a*x+s) - c*f(a*x-s)
-        # where we express a positive and a negative shift by default
-        default_param_shift = [
-            [multiplier, a, shift],
-            [-multiplier, a, -shift],  # pylint: disable=invalid-unary-operand-type
-        ]
-        return default_param_shift
+        # We no longer assume any default parameter-shift rule to apply.
+        raise OperatorPropertyUndefined(
+            f"The operation {self.name} does not have a parameter-shift recipe defined."
+            " This error might occur if previously the two-term shift rule was assumed"
+            " silently. In this case, consider adding it explicitly to the operation."
+        )
 
     @property
     def parameter_frequencies(self):
@@ -1291,18 +1295,25 @@ class Operation(Operator):
         if self.num_params == 1:
             # if the operator has a single parameter, we can query the
             # generator, and if defined, use its eigenvalues.
+            try:
+                gen = qml.generator(self, format="observable")
+            except GeneratorUndefinedError as e:
+                raise ParameterFrequenciesUndefinedError(
+                    f"Operation {self.name} does not have parameter frequencies defined."
+                ) from e
 
             with warnings.catch_warnings():
                 warnings.filterwarnings(
                     action="ignore", message=r".+ eigenvalues will be computed numerically\."
                 )
-                eigvals = qml.eigvals(qml.generator(self, format="observable"))
+                eigvals = qml.eigvals(gen)
 
             eigvals = tuple(np.round(eigvals, 8))
             return qml.gradients.eigvals_to_frequencies(eigvals)
 
-        raise OperatorPropertyUndefined(
-            f"Operation {self.name} does not have parameter frequencies."
+        raise ParameterFrequenciesUndefinedError(
+            f"Operation {self.name} does not have parameter frequencies defined, "
+            "and parameter frequencies can not be computed as no generator is defined."
         )
 
     @property
@@ -1391,12 +1402,8 @@ class Operation(Operator):
 
         # check the grad_recipe validity
         if self.grad_recipe is None:
-            # default recipe for every parameter
+            # Make sure grad_recipe is an iterable of correct length instead of None
             self.grad_recipe = [None] * self.num_params
-        else:
-            assert (
-                len(self.grad_recipe) == self.num_params
-            ), "Gradient recipe must have one entry for each parameter!"
 
 
 class Channel(Operation, abc.ABC):
diff --git a/pennylane/ops/cv.py b/pennylane/ops/cv.py
index 55eb663b0..d682f501c 100644
--- a/pennylane/ops/cv.py
+++ b/pennylane/ops/cv.py
@@ -46,6 +46,9 @@ from pennylane import math as qml_math
 from .identity import Identity  #  pylint: disable=unused-import
 
 
+_two_term_shift_rule = [[0.5, 1, np.pi / 2], [-0.5, 1, -np.pi / 2]]
+
+
 def _rotation(phi, bare=False):
     r"""Utility function, returns the Heisenberg transformation of a phase rotation gate.
 
@@ -104,6 +107,7 @@ class Rotation(CVOperation):
     num_params = 1
     num_wires = 1
     grad_method = "A"
+    grad_recipe = (_two_term_shift_rule,)
 
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
@@ -158,7 +162,7 @@ class Squeezing(CVOperation):
     shift = 0.1
     multiplier = 0.5 / math.sinh(shift)
     a = 1
-    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], None)
+    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], _two_term_shift_rule)
 
     def __init__(self, r, phi, wires, do_queue=True, id=None):
         super().__init__(r, phi, wires=wires, do_queue=do_queue, id=id)
@@ -215,7 +219,7 @@ class Displacement(CVOperation):
     shift = 0.1
     multiplier = 0.5 / shift
     a = 1
-    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], None)
+    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], _two_term_shift_rule)
 
     def __init__(self, a, phi, wires, do_queue=True, id=None):
         super().__init__(a, phi, wires=wires, do_queue=do_queue, id=id)
@@ -274,6 +278,7 @@ class Beamsplitter(CVOperation):
     num_params = 2
     num_wires = 2
     grad_method = "A"
+    grad_recipe = (_two_term_shift_rule, _two_term_shift_rule)
 
     def __init__(self, theta, phi, wires, do_queue=True, id=None):
         super().__init__(theta, phi, wires=wires, do_queue=do_queue, id=id)
@@ -342,7 +347,7 @@ class TwoModeSqueezing(CVOperation):
     shift = 0.1
     multiplier = 0.5 / math.sinh(shift)
     a = 1
-    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], None)
+    grad_recipe = ([[multiplier, a, shift], [-multiplier, a, -shift]], _two_term_shift_rule)
 
     def __init__(self, r, phi, wires, do_queue=True, id=None):
         super().__init__(r, phi, wires=wires, do_queue=do_queue, id=id)
diff --git a/pennylane/ops/qubit/parametric_ops.py b/pennylane/ops/qubit/parametric_ops.py
index 9ddad2dff..231d00dd7 100644
--- a/pennylane/ops/qubit/parametric_ops.py
+++ b/pennylane/ops/qubit/parametric_ops.py
@@ -29,7 +29,6 @@ from pennylane.ops.qubit.non_parametric_ops import PauliX, PauliY, PauliZ, Hadam
 from pennylane.utils import expand, pauli_eigs
 from pennylane.wires import Wires
 
-
 INV_SQRT2 = 1 / math.sqrt(2)
 
 
@@ -1164,14 +1163,6 @@ class PauliRot(Operation):
         return PauliRot(-self.parameters[0], self.parameters[1], wires=self.wires)
 
 
-# Four term gradient recipe for controlled rotations
-c1 = INV_SQRT2 * (np.sqrt(2) + 1) / 4
-c2 = INV_SQRT2 * (np.sqrt(2) - 1) / 4
-a = np.pi / 2
-b = 3 * np.pi / 2
-four_term_grad_recipe = ([[c1, 1, a], [-c1, 1, -a], [-c2, 1, b], [c2, 1, -b]],)
-
-
 class CRX(Operation):
     r"""
     The controlled-RX operator
@@ -1218,7 +1209,6 @@ class CRX(Operation):
 
     basis = "X"
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe
     parameter_frequencies = [(0.5, 1.0)]
 
     def generator(self):
@@ -1368,7 +1358,6 @@ class CRY(Operation):
 
     basis = "Y"
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe
     parameter_frequencies = [(0.5, 1.0)]
 
     def generator(self):
@@ -1511,7 +1500,6 @@ class CRZ(Operation):
 
     basis = "Z"
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe
     parameter_frequencies = [(0.5, 1.0)]
 
     def generator(self):
@@ -1675,7 +1663,6 @@ class CRot(Operation):
     """int: Number of trainable parameters that the operator depends on."""
 
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe * 3
     parameter_frequencies = [(0.5, 1.0), (0.5, 1.0), (0.5, 1.0)]
 
     def __init__(self, phi, theta, omega, wires, do_queue=True, id=None):
diff --git a/pennylane/ops/qubit/qchem_ops.py b/pennylane/ops/qubit/qchem_ops.py
index 882249d15..b5f74ac8c 100644
--- a/pennylane/ops/qubit/qchem_ops.py
+++ b/pennylane/ops/qubit/qchem_ops.py
@@ -16,22 +16,12 @@ This submodule contains the discrete-variable quantum operations that come
 from quantum chemistry applications.
 """
 # pylint:disable=abstract-method,arguments-differ,protected-access
-import math
 import numpy as np
 from scipy.sparse import coo_matrix
 
 import pennylane as qml
 from pennylane.operation import Operation
 
-INV_SQRT2 = 1 / math.sqrt(2)
-
-# Four term gradient recipe for controlled rotations
-c1 = INV_SQRT2 * (np.sqrt(2) + 1) / 4
-c2 = INV_SQRT2 * (np.sqrt(2) - 1) / 4
-a = np.pi / 2
-b = 3 * np.pi / 2
-four_term_grad_recipe = ([[c1, 1, a], [-c1, 1, -a], [-c2, 1, b], [c2, 1, -b]],)
-
 
 class SingleExcitation(Operation):
     r"""
@@ -90,10 +80,8 @@ class SingleExcitation(Operation):
     grad_method = "A"
     """Gradient computation method."""
 
-    grad_recipe = four_term_grad_recipe
-    """Gradient recipe for the parameter-shift method."""
-
     parameter_frequencies = [(0.5, 1.0)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         w1, w2 = self.wires
@@ -206,6 +194,7 @@ class SingleExcitationMinus(Operation):
     """Gradient computation method."""
 
     parameter_frequencies = [(1,)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         w1, w2 = self.wires
@@ -345,6 +334,7 @@ class SingleExcitationPlus(Operation):
     """Gradient computation method."""
 
     parameter_frequencies = [(1,)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         w1, w2 = self.wires
@@ -507,10 +497,8 @@ class DoubleExcitation(Operation):
     grad_method = "A"
     """Gradient computation method."""
 
-    grad_recipe = four_term_grad_recipe
-    """Gradient recipe for the parameter-shift method."""
-
     parameter_frequencies = [(0.5, 1.0)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         w0, w1, w2, w3 = self.wires
@@ -688,6 +676,7 @@ class DoubleExcitationPlus(Operation):
     """Gradient computation method."""
 
     parameter_frequencies = [(1,)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         G = -1 * np.eye(16, dtype=np.complex64)
@@ -783,6 +772,7 @@ class DoubleExcitationMinus(Operation):
     """Gradient computation method."""
 
     parameter_frequencies = [(1,)]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         G = np.eye(16, dtype=np.complex64)
@@ -900,24 +890,7 @@ class OrbitalRotation(Operation):
     """Gradient computation method."""
 
     parameter_frequencies = [(0.5, 1.0, 1.5, 2.0)]
-
-    @property
-    def grad_recipe(self):
-        r"""tuple(list[list[float]]): Gradient recipe for the
-        parameter-shift method.
-
-        This is a tuple with one nested list per operation parameter. For
-        parameter :math:`\phi_k`, the nested list contains elements of the form
-        :math:`[c_i, a_i, s_i]` where :math:`i` is the index of the
-        term, resulting in a gradient recipe of
-
-        .. math:: \frac{\partial}{\partial\phi_k}f = \sum_{i} c_i f(a_i \phi_k + s_i).
-
-        Since the ``OrbitalRotation`` operation has four parameter frequencies, this
-        corresponds to a parameter-shift rule with eight terms.
-        """
-        coeffs, shifts = qml.gradients.generate_shift_rule(self.parameter_frequencies[0])
-        return [np.stack([coeffs, np.ones_like(coeffs), shifts]).T]
+    """Frequencies of the operation parameter with respect to an expectation value."""
 
     def generator(self):
         w0, w1, w2, w3 = self.wires
diff --git a/pennylane/tape/cv_param_shift.py b/pennylane/tape/cv_param_shift.py
index 82949a0cf..eb5b4092d 100644
--- a/pennylane/tape/cv_param_shift.py
+++ b/pennylane/tape/cv_param_shift.py
@@ -27,7 +27,7 @@ import pennylane as qml
 from pennylane.measurements import MeasurementProcess
 from pennylane.tape import QuantumTape
 
-from .qubit_param_shift import QubitParamShiftTape
+from .qubit_param_shift import QubitParamShiftTape, _get_operation_recipe
 
 
 class CVParamShiftTape(QubitParamShiftTape):
@@ -230,16 +230,13 @@ class CVParamShiftTape(QubitParamShiftTape):
             tapes.
         """
 
-        t_idx = list(self.trainable_params)[idx]
-        op = self._par_info[t_idx]["op"]
-        p_idx = self._par_info[t_idx]["p_idx"]
-
-        param_shift = op.get_parameter_shift(p_idx)
+        op, p_idx = self.get_operation(idx)
+        param_shift = _get_operation_recipe(op, p_idx, None)
         shift = np.zeros_like(params)
 
         coeffs = []
         tapes = []
-        for c, _a, s in param_shift:
+        for c, _a, s in zip(*param_shift):
 
             shift[idx] = s
 
@@ -281,15 +278,12 @@ class CVParamShiftTape(QubitParamShiftTape):
             tapes.
         """
 
-        t_idx = list(self.trainable_params)[idx]
-        op = self._par_info[t_idx]["op"]
-        p_idx = self._par_info[t_idx]["p_idx"]
+        op, p_idx = self.get_operation(idx)
+        param_shift = _get_operation_recipe(op, p_idx, None)
 
         dev_wires = options["dev_wires"]
 
-        param_shift = op.get_parameter_shift(p_idx)
-
-        if len(param_shift) != 2:
+        if len(param_shift[0]) != 2:
             # The 2nd order CV parameter-shift rule only accepts two-term shifts
             raise NotImplementedError(
                 "Taking the analytic gradient for order-2 operators is "
@@ -297,8 +291,9 @@ class CVParamShiftTape(QubitParamShiftTape):
                 "gradient recipe of more than two terms."
             )
 
-        c1, a1, s1 = param_shift[0]
-        c2, a2, s2 = param_shift[1]
+        c1, c2 = param_shift[0]
+        a1, a2 = param_shift[1]
+        s1, s2 = param_shift[2]
 
         shift = np.zeros_like(params)
         shift[idx] = s1
diff --git a/pennylane/tape/qubit_param_shift.py b/pennylane/tape/qubit_param_shift.py
index f5cb9452b..8a4898904 100644
--- a/pennylane/tape/qubit_param_shift.py
+++ b/pennylane/tape/qubit_param_shift.py
@@ -29,6 +29,56 @@ from pennylane.math import toarray
 from .jacobian_tape import JacobianTape
 
 
+def _get_operation_recipe(op, p_idx, shifts):
+    """Utility function to return the parameter-shift rule
+    of the operation corresponding to trainable parameter
+    t_idx on tape. This is a copy from qml.gradients that
+    is replicated here to avoid circular dependencies, and
+    that can soon be removed together with QubitParamShiftTape.
+
+    This function performs multiple attempts to obtain the recipe:
+
+    - If the operation has a custom ``grad_recipe`` defined, it is used.
+
+    - If ``parameter_frequencies`` yield a result, the frequencies are
+      used to construct the general parameter-shift rule via
+      ``qml.gradients.generate_shift_rule``
+      Note that by default, the generator is used to compute the frequencies
+      if they are not provided by a custom implementation.
+
+    That is, the order of precedence is ``grad_recipe``, custom
+    ``parameter_frequencies`` and finally ``generator`` via the default
+    implementation of the frequencies.
+    """
+    # Try to use the stored grad_recipe of the operation
+    recipe = op.grad_recipe[p_idx]
+    if recipe is not None:
+        recipe = np.array(recipe).T
+        # remove all small coefficients and shifts
+        recipe[np.abs(recipe) < 1e-10] = 0
+
+        # remove columns where the coefficients are 0
+        recipe = recipe[:, ~(recipe[0] == 0)]
+        # sort columns according to abs(shift)
+        return recipe[:, np.argsort(np.abs(recipe)[-1])]
+
+    # Try to obtain frequencies, either via custom implementation or from generator eigvals
+    try:
+        frequencies = op.parameter_frequencies[p_idx]
+    except qml.operation.ParameterFrequenciesUndefinedError as e:
+        raise qml.operation.OperatorPropertyUndefined(
+            f"The operation {op.name} does not have a grad_recipe, parameter_frequencies or "
+            "a generator defined. No parameter shift rule can be applied."
+        ) from e
+
+    # Create shift rule from frequencies with given shifts
+    coeffs, shifts = qml.gradients.generate_shift_rule(frequencies, shifts=shifts, order=1)
+    # The generated shift rules do not include a rescaling of the parameter, only shifts.
+    mults = np.ones_like(coeffs)
+
+    return coeffs, mults, shifts
+
+
 class QubitParamShiftTape(JacobianTape):
     r"""Quantum tape for qubit parameter-shift analytic differentiation method.
 
@@ -146,18 +196,15 @@ class QubitParamShiftTape(JacobianTape):
             in addition to a post-processing function to be applied to the evaluated
             tapes.
         """
-        t_idx = list(self.trainable_params)[idx]
-        op = self._par_info[t_idx]["op"]
-        p_idx = self._par_info[t_idx]["p_idx"]
-
-        s = options.get("shift", np.pi / 2)
-        param_shift = op.get_parameter_shift(p_idx, shift=s)
+        # pylint: disable=unused-argument
+        op, p_idx = self.get_operation(idx)
+        param_shift = _get_operation_recipe(op, p_idx, None)
 
         shift = np.zeros_like(params)
         coeffs = []
         tapes = []
 
-        for c, a, s in param_shift:
+        for c, a, s in zip(*param_shift):
             shift[idx] = s
             shifted_tape = self.copy(copy_operations=True, tape_cls=QuantumTape)
             shifted_tape.set_parameters(a * params + shift)
@@ -343,16 +390,12 @@ class QubitParamShiftTape(JacobianTape):
         else:
             # No optimizations can be made, generate all 4 tapes
             for idx in idxs:
-                t_idx = list(self.trainable_params)[idx]
-                op = self._par_info[t_idx]["op"]
-                p_idx = self._par_info[t_idx]["p_idx"]
-                s = idx_shifts[idx]
-
-                param_shift = op.get_parameter_shift(p_idx, shift=s)
+                op, p_idx = self.get_operation(idx)
+                param_shift = _get_operation_recipe(op, p_idx, None)
                 param_shifts.append(param_shift)
 
-            for c1, a, s1 in param_shifts[0]:
-                for c2, _, s2 in param_shifts[1]:
+            for c1, a, s1 in zip(*param_shifts[0]):
+                for c2, _, s2 in zip(*param_shifts[1]):
                     c = c1 * c2
                     s = s1 * shift[i] + s2 * shift[j]
                     shifted_tape = self.copy(copy_operations=True, tape_cls=QuantumTape)
@@ -443,10 +486,9 @@ class QubitParamShiftTape(JacobianTape):
 
             # if this variable uses parameter-shift
             if grad_info["grad_method"] == "A":
-
-                # get_parameter_shift returns operation specific derivative formula
-                # for ops with 4-term gradients, the array will contain 4 pieces of data.
-                num_executions += len(grad_info["op"].get_parameter_shift(grad_info["p_idx"]))
+                op = grad_info["op"]
+                p_idx = grad_info["p_idx"]
+                num_executions += len(_get_operation_recipe(op, p_idx, None))
 
         info["num_parameter_shift_executions"] = num_executions
 
diff --git a/pennylane/templates/subroutines/commuting_evolution.py b/pennylane/templates/subroutines/commuting_evolution.py
index 749d086b7..cbc5017b6 100644
--- a/pennylane/templates/subroutines/commuting_evolution.py
+++ b/pennylane/templates/subroutines/commuting_evolution.py
@@ -90,7 +90,7 @@ class CommutingEvolution(Operation):
             coeffs = [1, -1]
             obs = [qml.PauliX(0) @ qml.PauliY(1), qml.PauliY(0) @ qml.PauliX(1)]
             hamiltonian = qml.Hamiltonian(coeffs, obs)
-            frequencies = (2,4)
+            frequencies = (2, 4)
 
             @qml.qnode(dev)
             def circuit(time):
diff --git a/pennylane/templates/subroutines/fermionic_double_excitation.py b/pennylane/templates/subroutines/fermionic_double_excitation.py
index 905b529ba..e8628fd26 100644
--- a/pennylane/templates/subroutines/fermionic_double_excitation.py
+++ b/pennylane/templates/subroutines/fermionic_double_excitation.py
@@ -15,23 +15,12 @@ r"""
 Contains the FermionicDoubleExcitation template.
 """
 # pylint: disable-msg=too-many-branches,too-many-arguments,protected-access
-import math
-
 import numpy as np
 import pennylane as qml
 from pennylane.operation import Operation, AnyWires
 from pennylane.ops import RZ, RX, CNOT, Hadamard
 
 
-# Four term gradient recipe for controlled rotations
-INV_SQRT2 = 1 / math.sqrt(2)
-c1 = INV_SQRT2 * (np.sqrt(2) + 1) / 4
-c2 = INV_SQRT2 * (np.sqrt(2) - 1) / 4
-a = np.pi / 2
-b = 3 * np.pi / 2
-four_term_grad_recipe = ([[c1, 1, a], [-c1, 1, -a], [-c2, 1, b], [c2, 1, -b]],)
-
-
 def _layer1(weight, s, r, q, p, set_cnot_wires):
     r"""Implement the first layer of the circuit to exponentiate the double-excitation
     operator entering the UCCSD ansatz.
@@ -505,7 +494,7 @@ class FermionicDoubleExcitation(Operation):
 
     num_wires = AnyWires
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe
+    parameter_frequencies = [(0.5, 1.0)]
 
     def __init__(self, weight, wires1=None, wires2=None, do_queue=True, id=None):
 
diff --git a/pennylane/templates/subroutines/fermionic_single_excitation.py b/pennylane/templates/subroutines/fermionic_single_excitation.py
index d185d0492..d8dd153ce 100644
--- a/pennylane/templates/subroutines/fermionic_single_excitation.py
+++ b/pennylane/templates/subroutines/fermionic_single_excitation.py
@@ -15,23 +15,12 @@ r"""
 Contains the FermionicSingleExcitation template.
 """
 # pylint: disable-msg=too-many-branches,too-many-arguments,protected-access
-import math
-
 import pennylane as qml
 from pennylane import numpy as np
 from pennylane.operation import Operation, AnyWires
 from pennylane.ops import RZ, RX, CNOT, Hadamard
 
 
-# Four term gradient recipe for controlled rotations
-INV_SQRT2 = 1 / math.sqrt(2)
-c1 = INV_SQRT2 * (np.sqrt(2) + 1) / 4
-c2 = INV_SQRT2 * (np.sqrt(2) - 1) / 4
-a = np.pi / 2
-b = 3 * np.pi / 2
-four_term_grad_recipe = ([[c1, 1, a], [-c1, 1, -a], [-c2, 1, b], [c2, 1, -b]],)
-
-
 class FermionicSingleExcitation(Operation):
     r"""Circuit to exponentiate the tensor product of Pauli matrices representing the
     single-excitation operator entering the Unitary Coupled-Cluster Singles
@@ -124,7 +113,7 @@ class FermionicSingleExcitation(Operation):
 
     num_wires = AnyWires
     grad_method = "A"
-    grad_recipe = four_term_grad_recipe
+    parameter_frequencies = [(0.5, 1.0)]
 
     def __init__(self, weight, wires=None, do_queue=True, id=None):
         if len(wires) < 2:


Test patch:
diff --git a/tests/gradients/test_parameter_shift.py b/tests/gradients/test_parameter_shift.py
index ef5d14187..6b70fa171 100644
--- a/tests/gradients/test_parameter_shift.py
+++ b/tests/gradients/test_parameter_shift.py
@@ -305,16 +305,18 @@ class TestParamShift:
         assert np.allclose(j1, [exp, 0])
         assert np.allclose(j2, [0, exp])
 
-    def test_grad_recipe_parameter_dependent(self):
+    def test_grad_recipe_parameter_dependent(self, monkeypatch):
         """Test that an operation with a gradient recipe that depends on
         its instantiated parameter values works correctly within the parameter
-        shift rule"""
+        shift rule. Also tests that grad_recipes supersedes paramter_frequencies.
+        """
 
-        class RX(qml.RX):
-            @property
-            def parameter_frequencies(self):
-                raise qml.operation.OperatorPropertyUndefined
+        def fail(*args, **kwargs):
+            raise qml.operation.ParameterFrequenciesUndefinedError
 
+        monkeypatch.setattr(qml.RX, "parameter_frequencies", fail)
+
+        class RX(qml.RX):
             @property
             def grad_recipe(self):
                 # The gradient is given by [f(2x) - f(0)] / (2 sin(x)), by subsituting
@@ -339,38 +341,38 @@ class TestParamShift:
         grad = fn(dev.batch_execute(tapes))
         assert np.allclose(grad, -np.sin(x))
 
-    @pytest.mark.parametrize("shifts", [None, [(0.1, 0.3, 0.5)]])
-    def test_default_two_term_rule(self, shifts):
-        """Test that the default two-term shift rule is used if no
-        other information is available to obtain the shift rule."""
+    def test_error_no_diff_info(self):
+        """Test that an error is raised if no grad_recipe, no parameter_frequencies
+        and no generator are found."""
 
         class RX(qml.RX):
+            """This copy of RX overwrites parameter_frequencies to report
+            missing information, disabling its differentiation."""
+
             @property
             def parameter_frequencies(self):
-                raise qml.operation.OperatorPropertyUndefined
+                raise qml.operation.ParameterFrequenciesUndefinedError
+
+        class NewOp(qml.operation.Operation):
+            """This new operation does not overwrite parameter_frequencies
+            but does not have a generator, disabling its differentiation."""
 
-            def generator(self):
-                raise GeneratorUndefinedError(f"Operation {self.name} does not have a generator")
+            num_params = 1
+            grad_method = "A"
+            num_wires = 1
 
         x = np.array(0.654, requires_grad=True)
         dev = qml.device("default.qubit", wires=2)
 
-        with qml.tape.JacobianTape() as tape:
-            RX(x, wires=0)
-            qml.expval(qml.PauliZ(0))
-
-        tapes, fn = qml.gradients.param_shift(tape, shifts=shifts)
-
-        assert len(tapes) == 2
-        if shifts is None:
-            assert tapes[0].operations[0].data[0] == x + np.pi / 2
-            assert tapes[1].operations[0].data[0] == x - np.pi / 2
-        else:
-            assert tapes[0].operations[0].data[0] == x + shifts[0][0]
-            assert tapes[1].operations[0].data[0] == x - shifts[0][0]
+        for op in [RX, NewOp]:
+            with qml.tape.JacobianTape() as tape:
+                op(x, wires=0)
+                qml.expval(qml.PauliZ(0))
 
-        grad = fn(dev.batch_execute(tapes))
-        assert np.allclose(grad, -np.sin(x))
+            with pytest.raises(
+                qml.operation.OperatorPropertyUndefined, match="does not have a grad_recipe"
+            ):
+                qml.gradients.param_shift(tape)
 
 
 class TestParameterShiftRule:
diff --git a/tests/tape/test_qnode_old.py b/tests/tape/test_qnode_old.py
index 06545288e..6f6e27e4d 100644
--- a/tests/tape/test_qnode_old.py
+++ b/tests/tape/test_qnode_old.py
@@ -151,14 +151,14 @@ class TestValidation:
             QNode._validate_backprop_method(dev, "another_interface")
 
     def test_parameter_shift_tape_qubit_device(self):
-        """Test that the get_parameter_shift_method method correctly and
+        """Test that the get_parameter_shift_tape method correctly and
         returns the correct tape for qubit devices."""
         dev = qml.device("default.qubit", wires=1)
         tape_class = QNode._get_parameter_shift_tape(dev)
         assert tape_class is QubitParamShiftTape
 
     def test_parameter_shift_tape_cv_device(self):
-        """Test that the get_parameter_shift_method method correctly and
+        """Test that the get_parameter_shift_tape method correctly and
         returns the correct tape for qubit devices."""
         dev = qml.device("default.gaussian", wires=1)
         tape_class = QNode._get_parameter_shift_tape(dev)
diff --git a/tests/tape/test_qubit_param_shift.py b/tests/tape/test_qubit_param_shift.py
index ef1b4289e..1a5a27b5c 100644
--- a/tests/tape/test_qubit_param_shift.py
+++ b/tests/tape/test_qubit_param_shift.py
@@ -18,6 +18,42 @@ import numpy as np
 import pennylane as qml
 from pennylane import numpy as pnp
 from pennylane.tape import QubitParamShiftTape
+from pennylane.tape.qubit_param_shift import _get_operation_recipe
+from pennylane.operation import (
+    Operation,
+    OperatorPropertyUndefined,
+    ParameterFrequenciesUndefinedError,
+)
+
+
+class TestGetOperationRecipe:
+    """Tests special cases for the _get_operation_recipe
+    copy in qubit_param_shift.py, the original is in
+    gradients/parameter_shift.py
+    """
+
+    class DummyOp0(Operation):
+        num_params = 1
+        num_wires = 1
+        grad_recipe = (None,)
+
+    class DummyOp1(Operation):
+        num_params = 1
+        num_wires = 1
+        grad_recipe = (None,)
+
+        @property
+        def parameter_frequencies(self):
+            raise ParameterFrequenciesUndefinedError
+
+    @pytest.mark.parametrize("Op", [DummyOp0, DummyOp1])
+    def test_error_no_grad_info(self, Op):
+        op = Op(0.1, wires=0)
+        with pytest.raises(
+            OperatorPropertyUndefined,
+            match=f"The operation {op.name} does not have a grad_recipe,",
+        ):
+            _get_operation_recipe(op, 0, None)
 
 
 class TestGradMethod:
diff --git a/tests/test_operation.py b/tests/test_operation.py
index ca3615a85..afaacb445 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -152,20 +152,6 @@ class TestOperatorConstruction:
 class TestOperationConstruction:
     """Test custom operations construction."""
 
-    def test_incorrect_grad_recipe_length(self):
-        """Test that an exception is raised if len(grad_recipe)!=len(num_params)"""
-
-        class DummyOp(qml.operation.CVOperation):
-            r"""Dummy custom operation"""
-            num_wires = 2
-            grad_method = "A"
-            grad_recipe = [(0.5, 0.1), (0.43, 0.1)]
-
-        with pytest.raises(
-            AssertionError, match="Gradient recipe must have one entry for each parameter"
-        ):
-            DummyOp(0.5, wires=[0, 1])
-
     def test_grad_recipe_parameter_dependent(self):
         """Test that an operation with a gradient recipe that depends on
         its instantiated parameter values works correctly"""
@@ -184,6 +170,37 @@ class TestOperationConstruction:
         op = DummyOp(x, wires=0)
         assert op.grad_recipe == ([[1.0, 1.0, x], [1.0, 0.0, -x]],)
 
+    def test_warning_get_parameter_shift(self):
+        """Test that ``get_parameter_shift`` issues a deprecation
+        warning."""
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+            num_params = 1
+            grad_recipe = ("Dummy recipe",)
+
+        op = DummyOp(0.1, wires=0)
+        with pytest.warns(UserWarning, match="get_parameter_shift is deprecated"):
+            assert op.get_parameter_shift(0) == "Dummy recipe"
+
+    def test_error_get_parameter_shift_no_recipe(self):
+        """Test that ``get_parameter_shift`` raises an Error if no grad_recipe
+        is available, as we no longer assume the two-term rule by default."""
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+            num_params = 1
+            grad_recipe = (None,)
+
+        op = DummyOp(0.1, wires=0)
+        with pytest.raises(
+            qml.operation.OperatorPropertyUndefined,
+            match="The operation DummyOp does not have a parameter-shift recipe",
+        ):
+            op.get_parameter_shift(0)
+
     def test_default_grad_method_with_frequencies(self):
         """Test that the correct ``grad_method`` is returned by default
         if ``parameter_frequencies`` are present.
@@ -230,6 +247,20 @@ class TestOperationConstruction:
         op = DummyOp(x, wires=0)
         assert op.grad_method == "F"
 
+    def test_default_grad_method_with_grad_recipe(self):
+        """Test that the correct ``grad_method`` is returned by default
+        if a grad_recipe is present.
+        """
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+            grad_recipe = ["not a recipe"]
+
+        x = 0.654
+        op = DummyOp(x, wires=0)
+        assert op.grad_method == "A"
+
     def test_default_grad_no_param(self):
         """Test that the correct ``grad_method`` is returned by default
         if an operation does not have a parameter.
@@ -271,7 +302,7 @@ class TestOperationConstruction:
         x = [0.654, 2.31, 0.1]
         op = DummyOp(*x, wires=0)
         with pytest.raises(
-            qml.operation.OperatorPropertyUndefined, match="does not have parameter"
+            qml.operation.OperatorPropertyUndefined, match="DummyOp does not have parameter"
         ):
             op.parameter_frequencies