From 1a6bb86ec3e4616bceb3bb19e8a61d996f00867e Mon Sep 17 00:00:00 2001 From: "Kevin J. Sung" Date: Thu, 21 Sep 2023 08:32:30 -0400 Subject: [PATCH] add num-num interaction and hopgate and improve gates docs --- python/ffsim/__init__.py | 2 + python/ffsim/gates/__init__.py | 2 + python/ffsim/gates/diag_coulomb.py | 3 +- python/ffsim/gates/gates.py | 151 ++++++++++++++- python/ffsim/gates/num_op_sum.py | 2 +- python/ffsim/gates/orbital_rotation.py | 2 +- tests/gates/gates_test.py | 247 +++++++++++++++++++++---- 7 files changed, 363 insertions(+), 46 deletions(-) diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py index 7f0abd67d..ceda0a504 100644 --- a/python/ffsim/__init__.py +++ b/python/ffsim/__init__.py @@ -33,7 +33,9 @@ from ffsim.gates import ( apply_diag_coulomb_evolution, apply_givens_rotation, + apply_hop_gate, apply_num_interaction, + apply_num_num_interaction, apply_num_op_prod_interaction, apply_num_op_sum_evolution, apply_orbital_rotation, diff --git a/python/ffsim/gates/__init__.py b/python/ffsim/gates/__init__.py index 60aee026b..51b886d80 100644 --- a/python/ffsim/gates/__init__.py +++ b/python/ffsim/gates/__init__.py @@ -13,7 +13,9 @@ from ffsim.gates.diag_coulomb import apply_diag_coulomb_evolution from ffsim.gates.gates import ( apply_givens_rotation, + apply_hop_gate, apply_num_interaction, + apply_num_num_interaction, apply_num_op_prod_interaction, apply_tunneling_interaction, ) diff --git a/python/ffsim/gates/diag_coulomb.py b/python/ffsim/gates/diag_coulomb.py index bfc682faa..2c4421c17 100644 --- a/python/ffsim/gates/diag_coulomb.py +++ b/python/ffsim/gates/diag_coulomb.py @@ -50,7 +50,8 @@ def apply_diag_coulomb_evolution( .. math:: \mathcal{U} - \exp(-i t \sum_{i, j, \sigma, \tau} Z_{ij} n_{i, \sigma} n_{j, \tau} / 2) + \exp\left(-i t \sum_{i, j, \sigma, \tau} + Z_{ij} n_{i, \sigma} n_{j, \tau} / 2\right) \mathcal{U}^\dagger where :math:`n_{i, \sigma}` denotes the number operator on orbital :math:`i` diff --git a/python/ffsim/gates/gates.py b/python/ffsim/gates/gates.py index d15d94b3d..5ec500f2a 100644 --- a/python/ffsim/gates/gates.py +++ b/python/ffsim/gates/gates.py @@ -62,7 +62,19 @@ def apply_givens_rotation( .. math:: - G(\theta) = \exp(\theta (a^\dagger_i a_j - a\dagger_j a_i)) + \text{G}(\theta) = \exp\left(\theta (a^\dagger_i a_j - a^\dagger_j a_i)\right) + + Under the Jordan-Wigner transform, this gate has the following matrix when applied + to neighboring qubits: + + .. math:: + + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + 0 & \cos(\theta) & -\sin(\theta) & 0\\ + 0 & \sin(\theta) & \cos(\theta) & 0\\ + 0 & 0 & 0 & 1 \\ + \end{pmatrix} Args: vec: The state vector to be transformed. @@ -102,7 +114,19 @@ def apply_tunneling_interaction( .. math:: - T(\theta) = \exp(i \theta (a^\dagger_i a_j + a\dagger_j a_i)) + \text{T}(\theta) = \exp\left(i \theta (a^\dagger_i a_j + a^\dagger_j a_i)\right) + + Under the Jordan-Wigner transform, this gate has the following matrix when applied + to neighboring qubits: + + .. math:: + + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + 0 & \cos(\theta) & i \sin(\theta) & 0\\ + 0 & i \sin(\theta) & \cos(\theta) & 0\\ + 0 & 0 & 0 & 1 \\ + \end{pmatrix} Args: vec: The state vector to be transformed. @@ -147,7 +171,7 @@ def apply_num_interaction( .. math:: - N(\theta) = \exp(i \theta a^\dagger_i a_i) + \text{N}(\theta) = \exp\left(i \theta a^\dagger_i a_i\right) Args: vec: The state vector to be transformed. @@ -184,6 +208,62 @@ def apply_num_interaction( return vec +def apply_num_num_interaction( + vec: np.ndarray, + theta: float, + target_orbs: tuple[int, int], + norb: int, + nelec: tuple[int, int], + *, + copy: bool = True, +): + r"""Apply a number-number interaction gate. + + The number-number interaction gate is + + .. math:: + + \text{NN}(\theta) = \exp\left(i \theta a^\dagger_i a_i a^\dagger_j a_j\right) + + Args: + vec: The state vector to be transformed. + theta: The rotation angle. + target_orbs: The orbitals on which to apply the interaction. + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + copy: Whether to copy the vector before operating on it. + - If ``copy=True`` then this function always returns a newly allocated + vector and the original vector is left untouched. + - If ``copy=False`` then this function may still return a newly allocated + vector, but the original vector may have its data overwritten. + It is also possible that the original vector is returned, + modified in-place. + """ + if len(set(target_orbs)) == 1: + raise ValueError( + f"The orbitals to interact must be distinct. Got {target_orbs}." + ) + if copy: + vec = vec.copy() + vec = apply_num_op_prod_interaction( + vec, + theta, + target_orbs=(target_orbs, []), + norb=norb, + nelec=nelec, + copy=False, + ) + vec = apply_num_op_prod_interaction( + vec, + theta, + target_orbs=([], target_orbs), + norb=norb, + nelec=nelec, + copy=False, + ) + return vec + + def apply_num_op_prod_interaction( vec: np.ndarray, theta: float, @@ -199,13 +279,15 @@ def apply_num_op_prod_interaction( .. math:: - NP(\theta) = \exp(i \theta \prod a^\dagger_{i, \sigma} a_{i, \sigma}) + \text{NP}(\theta) = + \exp\left(i \theta \prod a^\dagger_{i, \sigma} a_{i, \sigma}\right) Args: vec: The state vector to be transformed. theta: The rotation angle. - target_orbs: The orbitals on which to apply the interaction. This should - be a tuple of (orbital, spin) pairs. + target_orbs: A pair of lists of integers giving the orbitals on which to apply + the interaction. The first list specifies the alpha orbitals and the second + list specifies the beta orbitals. norb: The number of spatial orbitals. nelec: The number of alpha and beta electrons. copy: Whether to copy the vector before operating on it. @@ -228,3 +310,60 @@ def apply_num_op_prod_interaction( copy=False, ) return vec + + +def apply_hop_gate( + vec: np.ndarray, + theta: float, + target_orbs: tuple[int, int], + norb: int, + nelec: tuple[int, int], + *, + copy: bool = True, +) -> np.ndarray: + r"""Apply a hop gate. + + A "hop gate" is a Givens rotation gate followed by a number-number interaction + gate with angle pi: + + .. math:: + + \text{Hop}(\theta) = \text{NN}(\pi) \text{G}(\theta) + = \exp\left(i \pi a^\dagger_i a_i a^\dagger_j a_j\right) + \exp\left(\theta (a^\dagger_i a_j - a^\dagger_j a_i)\right) + + Under the Jordan-Wigner transform, this gate has the following matrix when applied + to neighboring qubits: + + .. math:: + + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + 0 & \cos(\theta) & -\sin(\theta) & 0\\ + 0 & \sin(\theta) & \cos(\theta) & 0\\ + 0 & 0 & 0 & -1 \\ + \end{pmatrix} + + Args: + vec: The state vector to be transformed. + theta: The rotation angle. + target_orbs: The orbitals (i, j) to rotate. + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + copy: Whether to copy the vector before operating on it. + - If ``copy=True`` then this function always returns a newly allocated + vector and the original vector is left untouched. + - If ``copy=False`` then this function may still return a newly allocated + vector, but the original vector may have its data overwritten. + It is also possible that the original vector is returned, + modified in-place. + """ + if copy: + vec = vec.copy() + vec = apply_givens_rotation( + vec, theta, target_orbs, norb=norb, nelec=nelec, copy=False + ) + vec = apply_num_num_interaction( + vec, np.pi, target_orbs, norb=norb, nelec=nelec, copy=False + ) + return vec diff --git a/python/ffsim/gates/num_op_sum.py b/python/ffsim/gates/num_op_sum.py index 30704a57f..0859683d6 100644 --- a/python/ffsim/gates/num_op_sum.py +++ b/python/ffsim/gates/num_op_sum.py @@ -41,7 +41,7 @@ def apply_num_op_sum_evolution( .. math:: \mathcal{U} - \exp(-i t \sum_{i, \sigma} \lambda_i n_{i, \sigma}) + \exp\left(-i t \sum_{i, \sigma} \lambda_i n_{i, \sigma}\right) \mathcal{U}^\dagger where :math:`n_{i, \sigma}` denotes the number operator on orbital :math:`i` diff --git a/python/ffsim/gates/orbital_rotation.py b/python/ffsim/gates/orbital_rotation.py index 035a790f8..15b6cc591 100644 --- a/python/ffsim/gates/orbital_rotation.py +++ b/python/ffsim/gates/orbital_rotation.py @@ -145,7 +145,7 @@ def apply_orbital_rotation( .. math:: - \exp(\sum_{ij} log(U)_{ij} a^\dagger{i} a_j) + \exp\left(\sum_{ij} \log(U)_{ij} a^\dagger_i a_j\right) Args: vec: The state vector to be transformed. diff --git a/tests/gates/gates_test.py b/tests/gates/gates_test.py index 74b4a1460..8b72cd838 100644 --- a/tests/gates/gates_test.py +++ b/tests/gates/gates_test.py @@ -13,6 +13,7 @@ from __future__ import annotations import itertools +from typing import Callable import numpy as np import pytest @@ -21,6 +22,56 @@ import ffsim +def assert_has_two_orbital_matrix( + gate: Callable[[np.ndarray, int, tuple[int, int]], np.ndarray], + target_orbs: tuple[int, int], + mat: np.ndarray, + phase_00: complex, + phase_11: complex, + norb: int, + rtol: float = 1e-7, + atol: float = 0.0, +): + state_00 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([], [])) + np.testing.assert_allclose( + np.vdot(state_00, gate(state_00, norb, (0, 0))), phase_00, rtol=rtol, atol=atol + ) + + i, j = target_orbs + + # test alpha + state_10 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([i], [])) + state_01 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([j], [])) + state_11 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([i, j], [])) + + np.testing.assert_allclose( + np.vdot(state_11, gate(state_11, norb, (2, 0))), phase_11, rtol=rtol, atol=atol + ) + + actual_mat = np.zeros((2, 2), dtype=complex) + for (a, state_a), (b, state_b) in itertools.product( + enumerate([state_01, state_10]), repeat=2 + ): + actual_mat[a, b] = np.vdot(state_a, gate(state_b, norb, (1, 0))) + np.testing.assert_allclose(actual_mat, mat, rtol=rtol, atol=atol) + + # test beta + state_10 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([], [i])) + state_01 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([], [j])) + state_11 = ffsim.slater_determinant(norb=norb, occupied_orbitals=([], [i, j])) + + np.testing.assert_allclose( + np.vdot(state_11, gate(state_11, norb, (0, 2))), phase_11, rtol=rtol, atol=atol + ) + + actual_mat = np.zeros((2, 2), dtype=complex) + for (a, state_a), (b, state_b) in itertools.product( + enumerate([state_01, state_10]), repeat=2 + ): + actual_mat[a, b] = np.vdot(state_a, gate(state_b, norb, (0, 1))) + np.testing.assert_allclose(actual_mat, mat, rtol=rtol, atol=atol) + + @pytest.mark.parametrize( "norb, nelec", [ @@ -30,29 +81,61 @@ ], ) def test_apply_givens_rotation(norb: int, nelec: tuple[int, int]): - """Test applying Givens rotation.""" + """Test Givens rotation.""" dim = ffsim.dim(norb, nelec) rng = np.random.default_rng() vec = np.array(ffsim.random.random_statevector(dim, seed=rng)) original_vec = vec.copy() - theta = rng.standard_normal() - for i, j in itertools.product(range(norb), repeat=2): - if i == j: - continue - result = ffsim.apply_givens_rotation(vec, theta, (i, j), norb=norb, nelec=nelec) - generator = np.zeros((norb, norb)) - generator[i, j] = theta - generator[j, i] = -theta - linop = ffsim.contract.hamiltonian_linop( - one_body_tensor=generator, norb=norb, nelec=nelec - ) - expected = scipy.sparse.linalg.expm_multiply( - linop, vec, traceA=np.sum(np.abs(generator)) - ) - np.testing.assert_allclose(result, expected, atol=1e-8) + theta = rng.uniform(-10, 10) + for i, j in itertools.combinations(range(norb), 2): + for target_orbs in [(i, j), (j, i)]: + result = ffsim.apply_givens_rotation( + vec, theta, target_orbs, norb=norb, nelec=nelec + ) + generator = np.zeros((norb, norb)) + a, b = target_orbs + generator[a, b] = theta + generator[b, a] = -theta + linop = ffsim.contract.hamiltonian_linop( + one_body_tensor=generator, norb=norb, nelec=nelec + ) + expected = scipy.sparse.linalg.expm_multiply( + linop, vec, traceA=np.sum(np.abs(generator)) + ) + np.testing.assert_allclose(result, expected, atol=1e-8) np.testing.assert_allclose(vec, original_vec) +def test_apply_givens_rotation_matrix(): + """Test Givens rotation matrix.""" + norb = 4 + rng = np.random.default_rng() + + def mat(theta: float) -> np.ndarray: + c = np.cos(theta) + s = np.sin(theta) + return np.array([[c, -s], [s, c]]) + + phase_00 = 1 + phase_11 = 1 + + for _ in range(5): + theta = 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, target_orbs=target_orbs, norb=norb, nelec=nelec + ), + target_orbs=target_orbs, + mat=mat(theta), + phase_00=phase_00, + phase_11=phase_11, + norb=norb, + atol=1e-12, + ) + + @pytest.mark.parametrize( "norb, nelec", [ @@ -62,27 +145,56 @@ def test_apply_givens_rotation(norb: int, nelec: tuple[int, int]): ], ) def test_apply_tunneling_interaction(norb: int, nelec: tuple[int, int]): - """Test applying tunneling interaction.""" + """Test tunneling interaction.""" dim = ffsim.dim(norb, nelec) rng = np.random.default_rng() vec = np.array(ffsim.random.random_statevector(dim, seed=rng)) - theta = rng.standard_normal() - for i, j in itertools.product(range(norb), repeat=2): - if i == j: - continue - result = ffsim.apply_tunneling_interaction( - vec, theta, (i, j), norb=norb, nelec=nelec - ) - generator = np.zeros((norb, norb)) - generator[i, j] = theta - generator[j, i] = theta - linop = ffsim.contract.hamiltonian_linop( - one_body_tensor=generator, norb=norb, nelec=nelec - ) - expected = scipy.sparse.linalg.expm_multiply( - 1j * linop, vec, traceA=np.sum(np.abs(generator)) - ) - np.testing.assert_allclose(result, expected, atol=1e-8) + theta = rng.uniform(-10, 10) + for i, j in itertools.combinations(range(norb), 2): + for target_orbs in [(i, j), (j, i)]: + result = ffsim.apply_tunneling_interaction( + vec, theta, target_orbs, norb=norb, nelec=nelec + ) + generator = np.zeros((norb, norb)) + generator[i, j] = theta + generator[j, i] = theta + linop = ffsim.contract.hamiltonian_linop( + one_body_tensor=generator, norb=norb, nelec=nelec + ) + expected = scipy.sparse.linalg.expm_multiply( + 1j * linop, vec, traceA=np.sum(np.abs(generator)) + ) + np.testing.assert_allclose(result, expected, atol=1e-8) + + +def test_apply_tunneling_interaction_matrix(): + """Test tunneling interaction matrix.""" + norb = 4 + rng = np.random.default_rng() + + def mat(theta: float) -> np.ndarray: + c = np.cos(theta) + s = np.sin(theta) + return np.array([[c, 1j * s], [1j * s, c]]) + + phase_00 = 1 + phase_11 = 1 + + for _ in range(5): + theta = 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_tunneling_interaction( + vec, theta, target_orbs=target_orbs, norb=norb, nelec=nelec + ), + target_orbs=target_orbs, + mat=mat(theta), + phase_00=phase_00, + phase_11=phase_11, + norb=norb, + atol=1e-12, + ) @pytest.mark.parametrize( @@ -98,7 +210,7 @@ def test_apply_num_interaction(norb: int, nelec: tuple[int, int]): dim = ffsim.dim(norb, nelec) rng = np.random.default_rng() vec = np.array(ffsim.random.random_statevector(dim, seed=rng)) - theta = rng.standard_normal() + theta = rng.uniform(-10, 10) for target_orb in range(norb): result = ffsim.apply_num_interaction( vec, theta, target_orb, norb=norb, nelec=nelec @@ -123,12 +235,44 @@ def test_apply_num_interaction(norb: int, nelec: tuple[int, int]): ], ) def test_apply_num_num_interaction(norb: int, nelec: tuple[int, int]): - """Test applying number-number interaction.""" + """Test applying number interaction.""" + rng = np.random.default_rng() + occupied_orbitals = ffsim.testing.random_occupied_orbitals(norb, nelec) + vec = ffsim.slater_determinant(norb, occupied_orbitals) + + theta = rng.uniform(-10, 10) + for i, j in itertools.combinations(range(norb), 2): + for target_orbs in [(i, j), (j, i)]: + result = ffsim.apply_num_num_interaction( + vec, + theta, + target_orbs=target_orbs, + norb=norb, + nelec=nelec, + ) + eig = 0.0 + for occ in occupied_orbitals: + if i in occ and j in occ: + eig += theta + expected = np.exp(1j * eig) * vec + np.testing.assert_allclose(result, expected, atol=1e-8) + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (2, (1, 1)), + (4, (2, 2)), + (5, (3, 2)), + ], +) +def test_apply_num_op_prod(norb: int, nelec: tuple[int, int]): + """Test applying number operator product interaction.""" rng = np.random.default_rng() occupied_orbitals = ffsim.testing.random_occupied_orbitals(norb, nelec) vec = ffsim.slater_determinant(norb, occupied_orbitals) - theta = rng.standard_normal() + theta = rng.uniform(-10, 10) for i, j in itertools.combinations_with_replacement(range(norb), 2): for spin_i, spin_j in itertools.product(range(2), repeat=2): target_orbs: tuple[set[int], set[int]] = (set(), set()) @@ -148,3 +292,32 @@ def test_apply_num_num_interaction(norb: int, nelec: tuple[int, int]): eig = 0 expected = np.exp(1j * eig) * vec np.testing.assert_allclose(result, expected, atol=1e-8) + + +def test_apply_hop_gate_matrix(): + """Test applying hop gate matrix.""" + norb = 4 + rng = np.random.default_rng() + + def mat(theta: float) -> np.ndarray: + c = np.cos(theta) + s = np.sin(theta) + return np.array([[c, -s], [s, c]]) + + phase_00 = 1 + phase_11 = -1 + + for _ in range(5): + theta = 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_hop_gate( + vec, theta, target_orbs=target_orbs, norb=norb, nelec=nelec + ), + target_orbs=target_orbs, + mat=mat(theta), + phase_00=phase_00, + phase_11=phase_11, + norb=norb, + )