diff --git a/qiskit_cold_atom/fermions/fermionic_state.py b/qiskit_cold_atom/fermions/fermionic_state.py index b2f5359..bfd15d8 100644 --- a/qiskit_cold_atom/fermions/fermionic_state.py +++ b/qiskit_cold_atom/fermions/fermionic_state.py @@ -41,7 +41,7 @@ def __init__(self, occupations: Union[List[int], List[List[int]]]): - If the occupations are not 0 or 1 """ - if isinstance(occupations[0], int): + if isinstance(occupations[0], (int, np.integer)): occupations = [occupations] self._sites = len(occupations[0]) @@ -153,7 +153,6 @@ def initial_state(cls, circuit: QuantumCircuit, num_species: int = 1) -> "Fermio # check that there are no more 'LoadFermions' instructions for instruction in circuit.data: - qargs = [circuit.qubits.index(qubit) for qubit in instruction[1]] if instruction[0].name == "load": diff --git a/qiskit_cold_atom/fermions/ffsim_backend.py b/qiskit_cold_atom/fermions/ffsim_backend.py index 01b3424..b98b49b 100644 --- a/qiskit_cold_atom/fermions/ffsim_backend.py +++ b/qiskit_cold_atom/fermions/ffsim_backend.py @@ -118,7 +118,6 @@ def _execute(self, data: Dict[str, Any], job_id: str = ""): output = {"results": []} num_species = data["num_species"] - assert num_species == 2 shots = data["shots"] seed = data["seed"] @@ -160,7 +159,7 @@ def _execute(self, data: Dict[str, Any], job_id: str = ""): if not measured_wires: shots = None - simulation_result = _simulate_ffsim(circuit, shots, seed) + simulation_result = _simulate_ffsim(circuit, num_species, shots, seed) output["results"].append( { @@ -189,7 +188,7 @@ def run( circuits: Union[QuantumCircuit, List[QuantumCircuit]], shots: int = 1000, seed: Optional[int] = None, - num_species: int = 2, + num_species: int = 1, **run_kwargs, ) -> AerJob: """ @@ -206,8 +205,6 @@ def run( Returns: aer_job: a job object containing the result of the simulation """ - assert num_species == 2 - if isinstance(circuits, QuantumCircuit): circuits = [circuits] @@ -229,10 +226,13 @@ def run( return aer_job -def _simulate_ffsim(circuit: QuantumCircuit, shots: int | None = None, seed=None) -> dict[str, Any]: +def _simulate_ffsim( + circuit: QuantumCircuit, num_species: int, shots: int | None = None, seed=None +) -> dict[str, Any]: assert circuit.num_qubits % 2 == 0 - norb = circuit.num_qubits // 2 - occ_a, occ_b = _get_initial_occupations(circuit) + assert num_species in (1, 2) + norb = circuit.num_qubits // num_species + occ_a, occ_b = _get_initial_occupations(circuit, num_species) nelec = len(occ_a), len(occ_b) vec = ffsim.slater_determinant(norb, (occ_a, occ_b)) qubit_indices = {q: i for i, q in enumerate(circuit.qubits)} @@ -240,26 +240,44 @@ def _simulate_ffsim(circuit: QuantumCircuit, shots: int | None = None, seed=None op, qubits, _ = instruction.operation, instruction.qubits, instruction.clbits if isinstance(op, Hop): orbs = [qubit_indices[q] for q in qubits] - spatial_orbs = _get_spatial_orbitals(orbs, norb) + spatial_orbs = _get_spatial_orbitals(orbs, norb, num_species) vec = _simulate_hop( - vec, np.array(op.params), spatial_orbs, norb=norb, nelec=nelec, copy=False + vec, + np.array(op.params), + spatial_orbs, + norb=norb, + nelec=nelec, + num_species=num_species, + copy=False, ) elif isinstance(op, Interaction): orbs = [qubit_indices[q] for q in qubits] - spatial_orbs = _get_spatial_orbitals(orbs, norb) + spatial_orbs = _get_spatial_orbitals(orbs, norb, num_species) (interaction,) = op.params vec = _simulate_interaction( - vec, interaction, spatial_orbs, norb=norb, nelec=nelec, copy=False + vec, + interaction, + spatial_orbs, + norb=norb, + nelec=nelec, + num_species=num_species, + copy=False, ) elif isinstance(op, Phase): orbs = [qubit_indices[q] for q in qubits] - spatial_orbs = _get_spatial_orbitals(orbs, norb) + spatial_orbs = _get_spatial_orbitals(orbs, norb, num_species) vec = _simulate_phase( - vec, np.array(op.params), spatial_orbs, norb=norb, nelec=nelec, copy=False + vec, + np.array(op.params), + spatial_orbs, + norb=norb, + nelec=nelec, + num_species=num_species, + copy=False, ) elif isinstance(op, FermionicGate): orbs = [qubit_indices[q] for q in qubits] - spatial_orbs = _get_spatial_orbitals(orbs, norb) + spatial_orbs = _get_spatial_orbitals(orbs, norb, num_species) ferm_op = _fermionic_op_to_fermion_operator(op.generator, spatial_orbs) linop = ffsim.linear_operator(ferm_op, norb, nelec) # TODO use ferm_op.values once it's available @@ -282,8 +300,9 @@ def _simulate_ffsim(circuit: QuantumCircuit, shots: int | None = None, seed=None return result -def _get_initial_occupations(circuit: QuantumCircuit): - norb = circuit.num_qubits // 2 +def _get_initial_occupations(circuit: QuantumCircuit, num_species: int): + assert num_species in (1, 2) + norb = circuit.num_qubits // num_species occ_a, occ_b = set(), set() occupations = [occ_a, occ_b] active_qubits = set() @@ -303,17 +322,63 @@ def _get_initial_occupations(circuit: QuantumCircuit): return tuple(occ_a), tuple(occ_b) -def _get_spatial_orbitals(orbs: list[int], norb: int) -> list[int]: - assert len(orbs) % 2 == 0 - alpha_orbs = orbs[: len(orbs) // 2] - beta_orbs = [orb - norb for orb in orbs[len(orbs) // 2 :]] - assert alpha_orbs == beta_orbs +def _get_spatial_orbitals(orbs: list[int], norb: int, num_species: int) -> list[int]: + assert num_species in (1, 2) + assert len(orbs) % num_species == 0 + alpha_orbs = orbs[: len(orbs) // num_species] + if num_species == 2: + beta_orbs = [orb - norb for orb in orbs[len(orbs) // 2 :]] + assert alpha_orbs == beta_orbs # reverse orbitals due to qiskit convention alpha_orbs = [norb - 1 - orb for orb in alpha_orbs] return alpha_orbs def _simulate_hop( + vec: np.ndarray, + coeffs: np.ndarray, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + num_species: int, + copy: bool, +) -> np.ndarray: + if num_species == 1: + return _simulate_hop_spinless( + vec=vec, coeffs=coeffs, target_orbs=target_orbs, norb=norb, nelec=nelec, copy=copy + ) + else: # num_species == 2 + return _simulate_hop_spinful( + vec=vec, coeffs=coeffs, target_orbs=target_orbs, norb=norb, nelec=nelec, copy=copy + ) + + +def _simulate_hop_spinless( + vec: np.ndarray, + coeffs: np.ndarray, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + copy: bool, +) -> np.ndarray: + assert norb % 2 == 0 + assert len(target_orbs) % 2 == 0 + mat = np.zeros((norb, norb)) + for i, val in zip(range(len(target_orbs) // 2 - 1), coeffs): + j, k = target_orbs[i], target_orbs[i + 1] + mat[j, k] = -val + mat[k, j] = -val + for i, val in zip(range(len(target_orbs) // 2, len(target_orbs) - 1), coeffs): + j, k = target_orbs[i], target_orbs[i + 1] + mat[j, k] = -val + mat[k, j] = -val + coeffs, orbital_rotation = scipy.linalg.eigh(mat) + return ffsim.apply_num_op_sum_evolution( + vec, coeffs, 1.0, norb=norb, nelec=nelec, orbital_rotation=orbital_rotation, copy=copy + ) + + +def _simulate_hop_spinful( vec: np.ndarray, coeffs: np.ndarray, target_orbs: list[int], @@ -333,6 +398,58 @@ def _simulate_hop( def _simulate_interaction( + vec: np.ndarray, + interaction: float, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + num_species: int, + copy: bool, +) -> np.ndarray: + if num_species == 1: + return _simulate_interaction_spinless( + vec=vec, + interaction=interaction, + target_orbs=target_orbs, + norb=norb, + nelec=nelec, + copy=copy, + ) + else: # num_species == 2 + return _simulate_interaction_spinful( + vec=vec, + interaction=interaction, + target_orbs=target_orbs, + norb=norb, + nelec=nelec, + copy=copy, + ) + + +def _simulate_interaction_spinless( + vec: np.ndarray, + interaction: float, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + copy: bool, +) -> np.ndarray: + assert len(target_orbs) % 2 == 0 + n_spatial_orbs = len(target_orbs) // 2 + mat = np.zeros((norb, norb)) + mat[target_orbs[:n_spatial_orbs], target_orbs[n_spatial_orbs:]] = interaction + mat[target_orbs[n_spatial_orbs:], target_orbs[:n_spatial_orbs]] = interaction + return ffsim.apply_diag_coulomb_evolution( + vec, + mat=mat, + time=1.0, + norb=norb, + nelec=nelec, + copy=copy, + ) + + +def _simulate_interaction_spinful( vec: np.ndarray, interaction: float, target_orbs: list[int], @@ -354,6 +471,53 @@ def _simulate_interaction( def _simulate_phase( + vec: np.ndarray, + mu: np.ndarray, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + num_species: int, + copy: bool, +) -> np.ndarray: + if num_species == 1: + return _simulate_phase_spinless( + vec=vec, + mu=mu, + target_orbs=target_orbs, + norb=norb, + nelec=nelec, + copy=copy, + ) + else: # num_species == 2 + return _simulate_phase_spinful( + vec=vec, + mu=mu, + target_orbs=target_orbs, + norb=norb, + nelec=nelec, + copy=copy, + ) + + +def _simulate_phase_spinless( + vec: np.ndarray, + mu: np.ndarray, + target_orbs: list[int], + norb: int, + nelec: tuple[int, int], + copy: bool, +) -> np.ndarray: + assert len(target_orbs) % 2 == 0 + n_spatial_orbs = len(target_orbs) // 2 + coeffs = np.zeros(norb) + coeffs[target_orbs[:n_spatial_orbs]] = mu + coeffs[target_orbs[n_spatial_orbs:]] = mu + return ffsim.apply_num_op_sum_evolution( + vec, coeffs, time=1.0, norb=norb, nelec=nelec, copy=copy + ) + + +def _simulate_phase_spinful( vec: np.ndarray, mu: np.ndarray, target_orbs: list[int], diff --git a/test/test_ffsim_backend.py b/test/test_ffsim_backend.py index 752bd94..cac74c3 100644 --- a/test/test_ffsim_backend.py +++ b/test/test_ffsim_backend.py @@ -39,9 +39,9 @@ def _random_occupations(norb: int, nelec: tuple[int, int], seed=None): rng = np.random.default_rng(seed) n_alpha, n_beta = nelec - alpha_bits = np.zeros(norb) + alpha_bits = np.zeros(norb, dtype=int) alpha_bits[:n_alpha] = 1 - beta_bits = np.zeros(norb) + beta_bits = np.zeros(norb, dtype=int) beta_bits[:n_beta] = 1 rng.shuffle(alpha_bits) rng.shuffle(beta_bits) @@ -63,36 +63,68 @@ def _fidelity(counts1: dict[str, int], counts2: dict[str, int]) -> float: class TestFfsimBackend(QiskitTestCase): """Test FfsimBackend.""" + def test_hop_gate_spinless(self): + """Test hop gate, spinless.""" + norb = 6 + nelec = 3 + + sim_backend = FermionSimulator() + ffsim_backend = FfsimBackend() + + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, (nelec, 0), seed=rng)[0] + hopping = rng.standard_normal(norb // 2 - 1) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Hop(norb, hopping), list(range(norb))) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + qubits = rng.choice(np.arange(norb), norb - 2, replace=False) + qc.append(Hop(norb - 2, hopping[: (norb - 2) // 2 - 1]), list(qubits)) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + def test_hop_gate(self): """Test hop gate.""" norb = 5 nelec = (3, 2) - rng = np.random.default_rng() - occupations = _random_occupations(norb, nelec, seed=rng) - hopping = rng.standard_normal(norb - 1) - sim_backend = FermionSimulator() ffsim_backend = FfsimBackend() - qc = sim_backend.initialize_circuit(occupations) - qc.append(Hop(2 * norb, hopping), list(range(2 * norb))) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) - - # test acting on subset of orbitals - qc = sim_backend.initialize_circuit(occupations) - orbs = rng.choice(np.arange(norb), norb - 1, replace=False) - qubits = np.concatenate([orbs, orbs + norb]) - qc.append(Hop(2 * (norb - 1), hopping[: norb - 2]), list(qubits)) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, nelec, seed=rng) + hopping = rng.standard_normal(norb - 1) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Hop(2 * norb, hopping), list(range(2 * norb))) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + orbs = rng.choice(np.arange(norb), norb - 1, replace=False) + qubits = np.concatenate([orbs, orbs + norb]) + qc.append(Hop(2 * (norb - 1), hopping[: norb - 2]), list(qubits)) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) def test_hop_gate_sign(self): """Test hop gate correctly handles fermionic sign.""" @@ -108,108 +140,212 @@ def test_hop_gate_sign(self): orbs = np.array([1, 0]) qubits = np.concatenate([orbs, orbs + norb]) qc.append(Hop(4, hopping), list(qubits)) - job = sim_backend.run(qc, shots=10, seed=1234, num_species=2) + job = sim_backend.run(qc, num_species=2) expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc, shots=10, seed=1234) + job = ffsim_backend.run(qc, num_species=2) ffsim_vec = job.result().get_statevector() np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + def test_interaction_gate_spinless(self): + """Test interaction gate, spinless.""" + norb = 6 + nelec = 3 + + sim_backend = FermionSimulator() + ffsim_backend = FfsimBackend() + + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, (nelec, 0), seed=rng)[0] + interaction = rng.standard_normal() + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Interaction(norb, interaction), list(range(norb))) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + qubits = rng.choice(np.arange(norb), norb - 2, replace=False) + qc.append(Interaction(norb - 2, interaction), list(qubits)) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + def test_interaction_gate(self): """Test interaction gate.""" norb = 5 nelec = (3, 2) + sim_backend = FermionSimulator() + ffsim_backend = FfsimBackend() + rng = np.random.default_rng() - occupations = _random_occupations(norb, nelec, seed=rng) - interaction = rng.standard_normal() + for _ in range(3): + occupations = _random_occupations(norb, nelec, seed=rng) + interaction = rng.standard_normal() + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Interaction(2 * norb, interaction), list(range(2 * norb))) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + orbs = rng.choice(np.arange(norb), norb - 1, replace=False) + qubits = np.concatenate([orbs, orbs + norb]) + qc.append(Interaction(2 * (norb - 1), interaction), list(qubits)) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + def test_phase_gate_spinless(self): + """Test phase gate.""" + norb = 6 + nelec = 3 sim_backend = FermionSimulator() ffsim_backend = FfsimBackend() - qc = sim_backend.initialize_circuit(occupations) - qc.append(Interaction(2 * norb, interaction), list(range(2 * norb))) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) - # test acting on subset of orbitals - qc = sim_backend.initialize_circuit(occupations) - orbs = rng.choice(np.arange(norb), norb - 1, replace=False) - qubits = np.concatenate([orbs, orbs + norb]) - qc.append(Interaction(2 * (norb - 1), interaction), list(qubits)) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, (nelec, 0), seed=rng)[0] + mu = rng.standard_normal(norb // 2) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Phase(norb, mu), list(range(norb))) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + qubits = rng.choice(np.arange(norb), norb - 2, replace=False) + qc.append(Phase(norb - 2, mu[: (norb - 2) // 2]), list(qubits)) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) def test_phase_gate(self): """Test phase gate.""" norb = 5 nelec = (3, 2) + sim_backend = FermionSimulator() + ffsim_backend = FfsimBackend() + rng = np.random.default_rng() - occupations = _random_occupations(norb, nelec, seed=rng) - mu = rng.standard_normal(norb) + for _ in range(3): + occupations = _random_occupations(norb, nelec, seed=rng) + mu = rng.standard_normal(norb) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(Phase(2 * norb, mu), list(range(2 * norb))) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + orbs = rng.choice(np.arange(norb), norb - 1, replace=False) + qubits = np.concatenate([orbs, orbs + norb]) + qc.append(Phase(2 * (norb - 1), mu[: norb - 1]), list(qubits)) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + def test_fermi_hubbard_gate_spinless(self): + """Test Fermi-Hubbard gate.""" + norb = 6 + nelec = 3 sim_backend = FermionSimulator() ffsim_backend = FfsimBackend() - qc = sim_backend.initialize_circuit(occupations) - qc.append(Phase(2 * norb, mu), list(range(2 * norb))) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) - - # test acting on subset of orbitals - qc = sim_backend.initialize_circuit(occupations) - orbs = rng.choice(np.arange(norb), norb - 1, replace=False) - qubits = np.concatenate([orbs, orbs + norb]) - qc.append(Phase(2 * (norb - 1), mu[: norb - 1]), list(qubits)) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, (nelec, 0), seed=rng)[0] + hopping = rng.standard_normal(norb // 2 - 1) + interaction = rng.standard_normal() + mu = rng.standard_normal(norb // 2) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(FermiHubbard(norb, hopping, interaction, mu), list(range(norb))) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + qubits = rng.choice(np.arange(norb), norb - 2, replace=False) + qc.append( + FermiHubbard( + norb - 2, hopping[: (norb - 2) // 2 - 1], interaction, mu[: (norb - 2) // 2] + ), + list(qubits), + ) + job = sim_backend.run(qc) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) def test_fermi_hubbard_gate(self): """Test Fermi-Hubbard gate.""" norb = 5 nelec = (3, 2) - rng = np.random.default_rng() - occupations = _random_occupations(norb, nelec, seed=rng) - hopping = rng.standard_normal(norb - 1) - interaction = rng.standard_normal() - mu = rng.standard_normal(norb) - sim_backend = FermionSimulator() ffsim_backend = FfsimBackend() - qc = sim_backend.initialize_circuit(occupations) - qc.append(FermiHubbard(2 * norb, hopping, interaction, mu), list(range(2 * norb))) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) - - # test acting on subset of orbitals - qc = sim_backend.initialize_circuit(occupations) - orbs = rng.choice(np.arange(norb), norb - 1, replace=False) - qubits = np.concatenate([orbs, orbs + norb]) - qc.append( - FermiHubbard(2 * (norb - 1), hopping[: norb - 2], interaction, mu[: norb - 1]), - list(qubits), - ) - job = sim_backend.run(qc, num_species=2) - expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) - ffsim_vec = job.result().get_statevector() - np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + rng = np.random.default_rng() + for _ in range(3): + occupations = _random_occupations(norb, nelec, seed=rng) + hopping = rng.standard_normal(norb - 1) + interaction = rng.standard_normal() + mu = rng.standard_normal(norb) + + qc = sim_backend.initialize_circuit(occupations) + qc.append(FermiHubbard(2 * norb, hopping, interaction, mu), list(range(2 * norb))) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) + + # test acting on subset of orbitals + qc = sim_backend.initialize_circuit(occupations) + orbs = rng.choice(np.arange(norb), norb - 1, replace=False) + qubits = np.concatenate([orbs, orbs + norb]) + qc.append( + FermiHubbard(2 * (norb - 1), hopping[: norb - 2], interaction, mu[: norb - 1]), + list(qubits), + ) + job = sim_backend.run(qc, num_species=2) + expected_vec = job.result().get_statevector() + job = ffsim_backend.run(qc, num_species=2) + ffsim_vec = job.result().get_statevector() + np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) def test_fermi_hubbard_gate_simple(self): """Test a simple Fermi-Hubbard gate.""" @@ -227,7 +363,7 @@ def test_fermi_hubbard_gate_simple(self): qc.append(FermiHubbard(2 * norb, hopping, interaction, mu), list(range(2 * norb))) job = sim_backend.run(qc, num_species=2) expected_vec = job.result().get_statevector() - job = ffsim_backend.run(qc) + job = ffsim_backend.run(qc, num_species=2) ffsim_vec = job.result().get_statevector() np.testing.assert_allclose(ffsim_vec, expected_vec, atol=1e-12) @@ -256,7 +392,7 @@ def test_simulate(self): expected_vec = result.get_statevector() expected_counts = result.get_counts() - job = ffsim_backend.run(qc, shots=10000, seed=1234) + job = ffsim_backend.run(qc, shots=10000, seed=1234, num_species=2) result = job.result() ffsim_vec = result.get_statevector() ffsim_counts = result.get_counts()