add complex phase to givens rotation (#236)

python/ffsim/gates/basic_gates.py

@@ -59,6 +59,7 @@ def apply_givens_rotation(
     nelec: int | tuple[int, int],
     spin: Spin = Spin.ALPHA_AND_BETA,
     *,
+    phi: float = 0.0,
     copy: bool = True,
 ) -> np.ndarray:
     r"""Apply a Givens rotation gate.
@@ -67,9 +68,11 @@ def apply_givens_rotation(
 
     .. math::
 
-        \text{G}(\theta, (p, q)) = \prod_{\sigma}
+        \text{G}(\theta, \varphi, (p, q)) = \prod_{\sigma}
+        \exp\left(i\varphi a^\dagger_{\sigma, p} a_{\sigma, p}\right)
         \exp\left(\theta (a^\dagger_{\sigma, p} a_{\sigma, q}
         - a^\dagger_{\sigma, q} a_{\sigma, p})\right)
+        \exp\left(-i\varphi a^\dagger_{\sigma, p} a_{\sigma, p}\right)
 
     Under the Jordan-Wigner transform, this gate has the following matrix when applied
     to neighboring qubits:
@@ -78,8 +81,8 @@ def apply_givens_rotation(
 
         \begin{pmatrix}
             1 & 0 & 0 & 0 \\
-            0 & \cos(\theta) & -\sin(\theta) & 0\\
-            0 & \sin(\theta) & \cos(\theta) & 0\\
+            0 & \cos(\theta) & -e^{-i\varphi}\sin(\theta) & 0\\
+            0 & e^{i\varphi}\sin(\theta) & \cos(\theta) & 0\\
             0 & 0 & 0 & 1 \\
         \end{pmatrix}
 
@@ -97,6 +100,7 @@ def apply_givens_rotation(
             - To act on only spin beta, pass :const:`ffsim.Spin.BETA`.
             - To act on both spin alpha and spin beta, pass
               :const:`ffsim.Spin.ALPHA_AND_BETA` (this is the default value).
+        phi: The optional phase angle.
         copy: Whether to copy the vector before operating on it.
 
             - If `copy=True` then this function always returns a newly allocated
@@ -109,9 +113,9 @@ def apply_givens_rotation(
     if len(set(target_orbs)) == 1:
         raise ValueError(f"The orbitals to rotate must be distinct. Got {target_orbs}.")
     c = math.cos(theta)
-    s = math.sin(theta)
-    mat = np.eye(norb)
-    mat[np.ix_(target_orbs, target_orbs)] = [[c, s], [-s, c]]
+    s = cmath.exp(1j * phi) * math.sin(theta)
+    mat = np.eye(norb, dtype=complex)
+    mat[np.ix_(target_orbs, target_orbs)] = [[c, s], [-s.conjugate(), c]]
     if isinstance(nelec, int):
         return apply_orbital_rotation(vec, mat, norb=norb, nelec=nelec, copy=copy)
     return apply_orbital_rotation(

tests/python/gates/basic_gates_test.py

@@ -121,29 +121,31 @@ def test_apply_givens_rotation_matrix_spinful(norb: int, spin: ffsim.Spin):
     """Test Givens rotation matrix."""
     rng = np.random.default_rng()
 
-    def mat(theta: float) -> np.ndarray:
+    def mat(theta: float, phi: float) -> np.ndarray:
         c = np.cos(theta)
-        s = np.sin(theta)
-        return np.array([[c, -s], [s, c]])
+        s = np.exp(1j * phi) * np.sin(theta)
+        return np.array([[c, -s.conjugate()], [s, c]])
 
     phase_00 = 1
     phase_11 = 1
 
     for _ in range(3):
         theta = rng.uniform(-10, 10)
+        phi = rng.uniform(-10, 10)
         for i, j in itertools.combinations(range(norb), 2):
             for target_orbs in [(i, j), (j, i)]:
                 assert_has_two_orbital_matrix(
                     lambda vec, norb, nelec: ffsim.apply_givens_rotation(
                         vec,
                         theta,
+                        phi=phi,
                         target_orbs=target_orbs,
                         norb=norb,
                         nelec=nelec,
                         spin=spin,
                     ),
                     target_orbs=target_orbs,
-                    mat=mat(theta),
+                    mat=mat(theta, phi),
                     phase_00=phase_00,
                     phase_11=phase_11,
                     norb=norb,
@@ -156,34 +158,64 @@ def test_apply_givens_rotation_matrix_spinless(norb: int):
     """Test Givens rotation matrix, spinless."""
     rng = np.random.default_rng()
 
-    def mat(theta: float) -> np.ndarray:
+    def mat(theta: float, phi: float) -> np.ndarray:
         c = np.cos(theta)
-        s = np.sin(theta)
-        return np.array([[c, -s], [s, c]])
+        s = np.exp(1j * phi) * np.sin(theta)
+        return np.array([[c, -s.conjugate()], [s, c]])
 
     phase_00 = 1
     phase_11 = 1
 
     for _ in range(3):
         theta = rng.uniform(-10, 10)
+        phi = rng.uniform(-10, 10)
         for i, j in itertools.combinations(range(norb), 2):
             for target_orbs in [(i, j), (j, i)]:
                 assert_has_two_orbital_matrix_spinless(
                     lambda vec, norb, nelec: ffsim.apply_givens_rotation(
                         vec,
                         theta,
+                        phi=phi,
                         target_orbs=target_orbs,
                         norb=norb,
                         nelec=nelec,
                     ),
                     target_orbs=target_orbs,
-                    mat=mat(theta),
+                    mat=mat(theta, phi),
                     phase_00=phase_00,
                     phase_11=phase_11,
                     norb=norb,
                 )
 
 
+def test_apply_givens_rotation_definition():
+    """Test definition of complex Givens in terms of real Givens and phases."""
+    norb = 5
+    nelec = (3, 2)
+    rng = np.random.default_rng()
+    theta = rng.uniform(-10, 10)
+    phi = rng.uniform(-10, 10)
+    vec = ffsim.random.random_statevector(ffsim.dim(norb, nelec), seed=rng)
+
+    # apply complex givens rotation
+    result = ffsim.apply_givens_rotation(
+        vec, theta, phi=phi, target_orbs=(1, 2), norb=norb, nelec=nelec
+    )
+
+    # get expected result using real givens rotation and phases
+    expected = ffsim.apply_num_interaction(
+        vec, -phi, target_orb=1, norb=norb, nelec=nelec
+    )
+    expected = ffsim.apply_givens_rotation(
+        expected, theta, target_orbs=(1, 2), norb=norb, nelec=nelec
+    )
+    expected = ffsim.apply_num_interaction(
+        expected, phi, target_orb=1, norb=norb, nelec=nelec
+    )
+
+    np.testing.assert_allclose(result, expected)
+
+
 @pytest.mark.parametrize("norb, spin", ffsim.testing.generate_norb_spin(range(4)))
 def test_apply_tunneling_interaction_matrix_spinful(norb: int, spin: ffsim.Spin):
     """Test tunneling interaction matrix."""