Repository: PennyLaneAI/pennylane
Pull request #3911 (2023-03-15 19:38:49)

Avoid adding zero-terms to the Hamiltonian in `rydberg_drive`
**Context**
This PR updates the `rydberg_drive` function to account for trivial detuning or amplitude

**Changes**
* Updated the `detuning` argument of `rydberg_drive` to be a keyword argument.
* Added checks to see if `detuning` is non-trivial to avoid adding trivial terms to the Hamiltonian.
* Added checks to see if `amplitude` is non-trivial to avoid adding trivial terms to the Hamiltonian.
* Raise `ValueError` if both `amplitude` and `detuning` are trivial, as that would result in all coefficients being zero and the Hamiltonian being trivial.


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.
mudit2812: Question for reviewers: Should the `rydberg_drive` docstring be updated to mention that if the function is used to define local drive, amplitude and phase must be zero? I'm confused about whether that belongs in the `rydberg_drive` docstring or somewhere in the Braket plugin docs.
lillian542: > Should the rydberg_drive docstring be updated to mention that if the function is used to define local drive, amplitude and phase must be zero?

I'd say no. I think we're saying hardware specifications go in the braket plugin docs, plus a PennyLane demo about accessing the hardware. Hardware-device specific details don't go in the general doc-string.
codecov[bot]: ## [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3911?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> :exclamation: No coverage uploaded for pull request base (`call_on_rydberg_hamiltonian@17850d0`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#section-missing-base-commit).
> The diff coverage is `n/a`.

```diff
@@                      Coverage Diff                       @@
##             call_on_rydberg_hamiltonian    #3911   +/-   ##
==============================================================
  Coverage                               ?   99.70%           
==============================================================
  Files                                  ?      342           
  Lines                                  ?    29952           
  Branches                               ?        0           
==============================================================
  Hits                                   ?    29863           
  Misses                                 ?       89           
  Partials                               ?        0           
```



Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

lillian542: I don't know how much we care about names, but "Local Drive for Rydberg Hamiltonian" doesn't seem like a very descriptive name at this point, since is addresses both avoiding the amplitude term and detuning term if they are zero, and is only related to local drive in-so-far as local drive on a Aquila is limited to the detuning term. 

Seems like a more general "Avoid adding zero-terms to the Hamiltonian in rydberg_drive" might be more fitting

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index c625d0cca..48ddc5df3 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -16,6 +16,7 @@
   * A new user-facing `rydberg_drive` function is added, which returns a `RydbergHamiltonian` containing
     the Hamiltonian of the interaction between a driving laser field and a group of atoms.
   [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
+  [(#3911)](https://github.com/PennyLaneAI/pennylane/pull/3911)
 
 * Added `Operation.__truediv__` dunder method to be able to divide operators.
   [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
diff --git a/pennylane/pulse/rydberg_hamiltonian.py b/pennylane/pulse/rydberg_hamiltonian.py
index b9e22226a..b185eb767 100644
--- a/pennylane/pulse/rydberg_hamiltonian.py
+++ b/pennylane/pulse/rydberg_hamiltonian.py
@@ -128,7 +128,7 @@ def rydberg_drive(amplitude, phase, detuning, wires):
 
     .. math::
         \frac{1}{2} \Omega(t) \sum_{i \in \text{wires}} (\cos(\phi)\sigma_i^x - \sin(\phi)\sigma_i^y) -
-        \frac{1}{2} \delta(t) \sum_{i \in \text{wires}} \sigma_i^z
+        \delta(t) \sum_{i \in \text{wires}} \sigma_i^z
 
     where :math:`\Omega`, :math:`\delta` and :math:`\phi` correspond to the rabi frequency, detuning
     and phase of the laser, :math:`i` correspond to the wire index, and :math:`\sigma^\alpha` for
@@ -165,7 +165,7 @@ def rydberg_drive(amplitude, phase, detuning, wires):
     .. code-block:: python
 
         atom_coordinates = [[0, 0], [0, 4], [4, 0], [4, 4]]
-        wires = [1, 2, 3, 4]
+        wires = [0, 1, 2, 3]
         H_i = qml.pulse.rydberg_interaction(atom_coordinates, wires)
 
         amplitude = lambda p, t: p * jnp.sin(jnp.pi * t)
@@ -179,7 +179,7 @@ def rydberg_drive(amplitude, phase, detuning, wires):
     ParametrizedHamiltonian: terms=3
 
     The first two terms of the drive Hamiltonian ``H_d`` correspond to the first sum (the sine and cosine terms),
-    describing drive between the ground and excited states. The third term corresponds to the shift term
+    describing drive between the ground and excited states. The third term corresponding to the shift term
     due to detuning from resonance. This drive term corresponds to a global drive that acts on all 4 wires of
     the device.
 
@@ -196,9 +196,9 @@ def rydberg_drive(amplitude, phase, detuning, wires):
 
     >>> params = [2.4]
     >>> circuit(params)
-    Array(0.97137696, dtype=float32)
+    Array(0.94301294, dtype=float64)
     >>> jax.grad(circuit)(params)
-    [Array(0.10493923, dtype=float32)]
+    [Array(0.59484969, dtype=float64)]
 
     We can also create a Hamiltonian with multiple local drives. The following circuit corresponds to the
     evolution where an additional local drive acting on wires ``[0, 1]`` is added to the Hamiltonian:
@@ -218,29 +218,44 @@ def rydberg_drive(amplitude, phase, detuning, wires):
             qml.evolve(H)(params, t=[0, 10])
             return qml.expval(qml.PauliZ(0))
 
-    >>> params = [2.4, [1.3, -2.0]]
+    >>> params = [2.4, [1.3, -2.0], -1.5]
     >>> circuit_local(params)
-    Array(0.45782223, dtype=float64)
+    Array(0.80238028, dtype=float64)
     >>> jax.grad(circuit_local)(params)
-    [Array(-0.33522988, dtype=float64),
-     [Array(0.40320718, dtype=float64, weak_type=True),
-      Array(-0.12003976, dtype=float64, weak_type=True)]]
+    [Array(-0.73273312, dtype=float64),
+     [Array(-0.45635171, dtype=float64, weak_type=True),
+      Array(0.76572817, dtype=float64, weak_type=True)],
+     Array(0.01027312, dtype=float64)]
     """
     if isinstance(wires, int):
         wires = [wires]
+    trivial_detuning = not callable(detuning) and qml.math.isclose(detuning, 0.0)
 
-    # We compute the `coeffs` and `observables` of the laser field
-    coeffs = [
-        amplitude_and_phase(qml.math.cos, amplitude, phase),
-        amplitude_and_phase(qml.math.sin, amplitude, phase),
-        detuning,
-    ]
+    if not callable(amplitude) and qml.math.isclose(amplitude, 0.0):
+        if trivial_detuning:
+            raise ValueError(
+                f"Expected non-zero value for at least one of either amplitude or detuning, but received amplitude={amplitude} and detuning={detuning}. All terms are zero."
+            )
+
+        coeffs = []
+        observables = []
+
+    else:
+        # We compute the `coeffs` and `observables` of the laser field
+        coeffs = [
+            amplitude_and_phase(qml.math.cos, amplitude, phase),
+            amplitude_and_phase(qml.math.sin, amplitude, phase),
+        ]
+
+        drive_x_term = 0.5 * sum(qml.PauliX(wire) for wire in wires)
+        drive_y_term = -0.5 * sum(qml.PauliY(wire) for wire in wires)
 
-    drive_x_term = sum(qml.PauliX(wire) for wire in wires)
-    drive_y_term = sum(-qml.PauliY(wire) for wire in wires)
-    detuning_term = sum(qml.PauliZ(wire) for wire in wires)
+        observables = [drive_x_term, drive_y_term]
 
-    observables = [drive_x_term, drive_y_term, detuning_term]
+    if not trivial_detuning:
+        detuning_term = -1.0 * sum(qml.PauliZ(wire) for wire in wires)
+        coeffs.append(detuning)
+        observables.append(detuning_term)
 
     # We convert the pulse data into a list of ``RydbergPulse`` objects
     pulses = [RydbergPulse(amplitude, phase, detuning, wires)]


Test patch:
diff --git a/tests/pulse/test_rydberg_hamiltonian.py b/tests/pulse/test_rydberg_hamiltonian.py
index f687c29aa..45df34662 100644
--- a/tests/pulse/test_rydberg_hamiltonian.py
+++ b/tests/pulse/test_rydberg_hamiltonian.py
@@ -82,7 +82,10 @@ class TestRydbergHamiltonian:
             for op1, op2 in zip(sum_rm.ops, [qml.PauliX(4), qml.PauliZ(8), qml.PauliY(8)])
         )
         assert qml.math.allequal(sum_rm.register, atom_coordinates)
-        assert sum_rm.pulses == [RydbergPulse(1, 2, 3, [4, 8]), RydbergPulse(5, 6, 7, 8)]
+        assert sum_rm.pulses == [
+            RydbergPulse(1, 2, 3, [4, 8]),
+            RydbergPulse(5, 6, 7, 8),
+        ]
 
     def test_add_parametrized_hamiltonian(self):
         """Tests that adding a `RydbergHamiltonian` and `ParametrizedHamiltonian` works as
@@ -356,12 +359,12 @@ class TestRydbergDrive:
         Hd = H1 + H2
 
         ops_expected = [
-            qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(3)]),
-            qml.Hamiltonian([1, 1], [qml.PauliX(1), qml.PauliX(2)]),
-            qml.sum(qml.s_prod(-1, qml.PauliY(1)), qml.s_prod(-1, qml.PauliY(2))),
-            qml.Hamiltonian([1, 1], [qml.PauliX(0), qml.PauliX(3)]),
-            qml.sum(qml.s_prod(-1, qml.PauliY(0)), qml.s_prod(-1, qml.PauliY(3))),
-            qml.Hamiltonian([1, 1], [qml.PauliZ(1), qml.PauliZ(2)]),
+            qml.Hamiltonian([-1, -1], [qml.PauliZ(0), qml.PauliZ(3)]),
+            qml.Hamiltonian([0.5, 0.5], [qml.PauliX(1), qml.PauliX(2)]),
+            qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(1), qml.PauliY(2)]),
+            qml.Hamiltonian([0.5, 0.5], [qml.PauliX(0), qml.PauliX(3)]),
+            qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(0), qml.PauliY(3)]),
+            qml.Hamiltonian([-1, -1], [qml.PauliZ(1), qml.PauliZ(2)]),
         ]
         coeffs_expected = [
             3,
@@ -394,6 +397,60 @@ class TestRydbergDrive:
         # Hamiltonian is as expected
         assert qml.equal(Hd([0.5, -0.5], t=5), H_expected([0.5, -0.5], t=5))
 
+    def test_no_amplitude(self):
+        """Test that when amplitude is not specified, the drive term is correctly defined."""
+
+        def f(p, t):
+            return np.cos(p * t)
+
+        Hd = rydberg_drive(amplitude=0, phase=1, detuning=f, wires=[0, 3])
+
+        ops_expected = [qml.Hamiltonian([-1, -1], [qml.PauliZ(0), qml.PauliZ(3)])]
+        coeffs_expected = [f]
+        H_expected = RydbergHamiltonian(coeffs_expected, ops_expected)
+
+        assert qml.equal(Hd([0.1], 10), H_expected([0.1], 10))
+        assert isinstance(Hd, RydbergHamiltonian)
+        assert Hd.interaction_coeff == 862690
+        assert Hd.wires == Wires([0, 3])
+        assert Hd.register is None
+        assert len(Hd.coeffs) == 1
+        assert Hd.coeffs[0] is f
+        assert len(Hd.ops) == 1
+        assert qml.equal(Hd.ops[0], ops_expected[0])
+
+    def test_no_detuning(self):
+        """Test that when detuning not specified, the drive term is correctly defined."""
+
+        def f(p, t):
+            return np.cos(p * t)
+
+        Hd = rydberg_drive(amplitude=f, phase=1, detuning=0, wires=[0, 3])
+
+        ops_expected = [
+            qml.Hamiltonian([0.5, 0.5], [qml.PauliX(0), qml.PauliX(3)]),
+            qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(0), qml.PauliY(3)]),
+        ]
+        coeffs_expected = [
+            AmplitudeAndPhase(np.cos, f, 1),
+            AmplitudeAndPhase(np.sin, f, 1),
+        ]
+        H_expected = RydbergHamiltonian(coeffs_expected, ops_expected)
+
+        assert qml.equal(Hd([0.1], 10), H_expected([0.1], 10))
+        assert isinstance(Hd, RydbergHamiltonian)
+        assert Hd.interaction_coeff == 862690
+        assert Hd.wires == Wires([0, 3])
+        assert Hd.register is None
+        assert all(isinstance(coeff, AmplitudeAndPhase) for coeff in Hd.coeffs)
+        assert len(Hd.coeffs) == 2
+        assert all(qml.equal(op, op_expected) for op, op_expected in zip(Hd.ops, ops_expected))
+
+    def test_no_amplitude_no_detuning(self):
+        """Test that the correct error is raised if both amplitude and detuning are trivial."""
+        with pytest.raises(ValueError, match="Expected non-zero value for at least one of either"):
+            _ = rydberg_drive(0, np.pi, 0, wires=[0])
+
 
 def callable_amp(p, t):
     return np.polyval(p, t)
@@ -489,14 +546,14 @@ class TestAmplitudeAndPhase:
         evaluated_H = Hd([3.4, 5.6], t)
 
         expected_H_fixed = qml.s_prod(
-            detuning, qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(1)])
+            detuning, qml.Hamiltonian([-1, -1], [qml.PauliZ(0), qml.PauliZ(1)])
         )
 
         c1 = np.sin(3.4 * t) * np.cos(np.cos(5.6 * t))
         c2 = np.sin(3.4 * t) * np.sin(np.cos(5.6 * t))
         expected_H_parametrized = qml.sum(
-            qml.s_prod(c1, qml.Hamiltonian([1, 1], [qml.PauliX(0), qml.PauliX(1)])),
-            qml.s_prod(c2, qml.sum(qml.s_prod(-1, qml.PauliY(0)), qml.s_prod(-1, qml.PauliY(1)))),
+            qml.s_prod(c1, qml.Hamiltonian([0.5, 0.5], [qml.PauliX(0), qml.PauliX(1)])),
+            qml.s_prod(c2, qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(0), qml.PauliY(1)])),
         )
 
         assert qml.equal(evaluated_H[0], expected_H_fixed)
@@ -518,14 +575,14 @@ class TestAmplitudeAndPhase:
         evaluated_H = Hd([5.6], t)
 
         expected_H_fixed = qml.s_prod(
-            detuning, qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(1)])
+            detuning, qml.Hamiltonian([-1, -1], [qml.PauliZ(0), qml.PauliZ(1)])
         )
 
         c1 = 7.2 * np.cos(np.sin(5.6 * t))
         c2 = 7.2 * np.sin(np.sin(5.6 * t))
         expected_H_parametrized = qml.sum(
-            qml.s_prod(c1, qml.Hamiltonian([1, 1], [qml.PauliX(0), qml.PauliX(1)])),
-            qml.s_prod(c2, qml.sum(qml.s_prod(-1, qml.PauliY(0)), qml.s_prod(-1, qml.PauliY(1)))),
+            qml.s_prod(c1, qml.Hamiltonian([0.5, 0.5], [qml.PauliX(0), qml.PauliX(1)])),
+            qml.s_prod(c2, qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(0), qml.PauliY(1)])),
         )
 
         assert qml.equal(evaluated_H[0], expected_H_fixed)
@@ -547,14 +604,14 @@ class TestAmplitudeAndPhase:
         evaluated_H = Hd([3.4], t)
 
         expected_H_fixed = qml.s_prod(
-            detuning, qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(1)])
+            detuning, qml.Hamiltonian([-1, -1], [qml.PauliZ(0), qml.PauliZ(1)])
         )
 
         c1 = np.sin(3.4 * t) * np.cos(4.3)
         c2 = np.sin(3.4 * t) * np.sin(4.3)
         expected_H_parametrized = qml.sum(
-            qml.s_prod(c1, qml.Hamiltonian([1, 1], [qml.PauliX(0), qml.PauliX(1)])),
-            qml.s_prod(c2, qml.sum(qml.s_prod(-1, qml.PauliY(0)), qml.s_prod(-1, qml.PauliY(1)))),
+            qml.s_prod(c1, qml.Hamiltonian([0.5, 0.5], [qml.PauliX(0), qml.PauliX(1)])),
+            qml.s_prod(c2, qml.Hamiltonian([-0.5, -0.5], [qml.PauliY(0), qml.PauliY(1)])),
         )
 
         assert qml.equal(evaluated_H[0], expected_H_fixed)
@@ -679,8 +736,12 @@ class TestIntegration:
         def fc(p, t):
             return p[0] * jnp.sin(t) + jnp.cos(p[1] * t)
 
-        H1 = rydberg_drive(amplitude=fa, phase=0, detuning=fb, wires=1)
-        H2 = rydberg_drive(amplitude=fc, phase=3 * jnp.pi, detuning=jnp.pi / 4, wires=4)
+        def fd(p, t):
+            return p * jnp.cos(t)
+
+        H1 = rydberg_drive(amplitude=fa, phase=0, detuning=fb, wires=wires)
+        H2 = rydberg_drive(amplitude=fc, phase=3 * jnp.pi, detuning=0, wires=4)
+        H3 = rydberg_drive(amplitude=0, phase=0, detuning=fd, wires=[3, 0])
 
         dev = qml.device("default.qubit", wires=wires)
 
@@ -689,16 +750,21 @@ class TestIntegration:
 
         @qml.qnode(dev, interface="jax")
         def qnode(params):
-            qml.evolve(Hd + H1 + H2)(params, ts)
+            qml.evolve(Hd + H1 + H2 + H3)(params, ts)
             return qml.expval(H_obj)
 
         @jax.jit
         @qml.qnode(dev, interface="jax")
         def qnode_jit(params):
-            qml.evolve(Hd + H1 + H2)(params, ts)
+            qml.evolve(Hd + H1 + H2 + H3)(params, ts)
             return qml.expval(H_obj)
 
-        params = (jnp.ones(5), jnp.array([1.0, jnp.pi]), jnp.array([jnp.pi / 2, 0.5]))
+        params = (
+            jnp.ones(5),
+            jnp.array([1.0, jnp.pi]),
+            jnp.array([jnp.pi / 2, 0.5]),
+            jnp.array(-0.5),
+        )
         res = qnode(params)
         res_jit = qnode_jit(params)
 
diff --git a/tests/returntypes/qnn/test_keras_new.py b/tests/returntypes/qnn/test_keras_new.py
index db3bdad29..89ad27e4c 100644
--- a/tests/returntypes/qnn/test_keras_new.py
+++ b/tests/returntypes/qnn/test_keras_new.py
@@ -23,6 +23,8 @@ KerasLayer = qml.qnn.keras.KerasLayer
 
 tf = pytest.importorskip("tensorflow", minversion="2")
 
+# pylint: disable=unnecessary-dunder-call
+
 
 @pytest.fixture
 def model(get_circuit, n_qubits, output_dim):
diff --git a/tests/returntypes/qnn/test_qnn_torch_new.py b/tests/returntypes/qnn/test_qnn_torch_new.py
index 2f52a95cf..a7e8e61cd 100644
--- a/tests/returntypes/qnn/test_qnn_torch_new.py
+++ b/tests/returntypes/qnn/test_qnn_torch_new.py
@@ -24,6 +24,8 @@ from pennylane.qnn.torch import TorchLayer
 
 torch = pytest.importorskip("torch")
 
+# pylint: disable=unnecessary-dunder-call
+
 
 def indices_up_to(n_max):
     """Returns an iterator over the number of qubits and output dimension, up to value n_max.