Source code for ffsim.qiskit.gates.diag_coulomb_trotter

# (C) Copyright IBM 2024.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Diagonal Coulomb Trotter evolution gate."""

from __future__ import annotations

from collections.abc import Generator, Iterator, Sequence

import numpy as np
import scipy.linalg
from qiskit.circuit import (
    CircuitInstruction,
    Gate,
    QuantumCircuit,
    QuantumRegister,
    Qubit,
)
from qiskit.circuit.library import GlobalPhaseGate

from ffsim.hamiltonians import DiagonalCoulombHamiltonian
from ffsim.qiskit.gates.diag_coulomb import DiagCoulombEvolutionJW
from ffsim.qiskit.gates.num_op_sum import NumOpSumEvolutionJW
from ffsim.qiskit.gates.orbital_rotation import OrbitalRotationJW
from ffsim.trotter._util import simulate_trotter_step_iterator


[docs] class SimulateTrotterDiagCoulombSplitOpJW(Gate): r"""Split operator Trotter evolution of diagonal Coulomb Hamiltonian, Jordan-Wigner. This gate assumes that qubits are ordered such that the first `norb` qubits correspond to the alpha orbitals and the last `norb` qubits correspond to the beta orbitals. """
[docs] def __init__( self, hamiltonian: DiagonalCoulombHamiltonian, time: float, *, n_steps: int = 1, order: int = 0, label: str | None = None, ): r"""Create diagonal Coulomb split-operator Trotter evolution gate. Args: norb: The number of spatial orbitals. hamiltonian: The Hamiltonian. time: The evolution time. n_steps: The number of Trotter steps. order: The order of the Trotter decomposition. label: The label of the gate. """ if order < 0: raise ValueError(f"order must be non-negative, got {order}.") if n_steps < 0: raise ValueError(f"n_steps must be non-negative, got {n_steps}.") self.hamiltonian = hamiltonian self.time = time self.n_steps = n_steps self.order = order super().__init__( "dc_trotter_split_op_jw", 2 * self.hamiltonian.norb, [], label=label )
def _define(self): """Gate decomposition.""" qubits = QuantumRegister(self.num_qubits) self.definition = QuantumCircuit.from_instructions( _simulate_trotter_diag_coulomb_split_op( qubits, self.hamiltonian, time=self.time, n_steps=self.n_steps, order=self.order, ), qubits=qubits, )
def _simulate_trotter_diag_coulomb_split_op( qubits: Sequence[Qubit], hamiltonian: DiagonalCoulombHamiltonian, time: float, n_steps: int = 1, order: int = 0, ) -> Iterator[CircuitInstruction]: if n_steps == 0: return one_body_energies, one_body_basis_change = scipy.linalg.eigh( hamiltonian.one_body_tensor ) step_time = time / n_steps current_basis = np.eye(hamiltonian.norb, dtype=complex) for _ in range(n_steps): current_basis = yield from _simulate_trotter_step_diag_coulomb_split_op( qubits, current_basis, one_body_energies, one_body_basis_change, hamiltonian.diag_coulomb_mats, step_time, norb=hamiltonian.norb, order=order, ) yield CircuitInstruction(OrbitalRotationJW(hamiltonian.norb, current_basis), qubits) yield CircuitInstruction(GlobalPhaseGate(-time * hamiltonian.constant), []) def _simulate_trotter_step_diag_coulomb_split_op( qubits: Sequence[Qubit], current_basis: np.ndarray, one_body_energies: np.ndarray, one_body_basis_change: np.ndarray, diag_coulomb_mats: np.ndarray, time: float, norb: int, order: int, ) -> Generator[CircuitInstruction, None, np.ndarray]: diag_coulomb_aa, diag_coulomb_ab = diag_coulomb_mats eye = np.eye(norb) for term_index, time in simulate_trotter_step_iterator(2, time, order): if term_index == 0: yield CircuitInstruction( OrbitalRotationJW(norb, one_body_basis_change.T.conj() @ current_basis), qubits, ) yield CircuitInstruction( NumOpSumEvolutionJW(norb, coeffs=one_body_energies, time=time), qubits ) current_basis = one_body_basis_change else: yield CircuitInstruction(OrbitalRotationJW(norb, current_basis), qubits) yield CircuitInstruction( DiagCoulombEvolutionJW( norb, (diag_coulomb_aa, diag_coulomb_ab, diag_coulomb_aa), time, ), qubits, ) current_basis = eye return current_basis