From 97646a3e716495e463c4baa46aca8399563a1837 Mon Sep 17 00:00:00 2001 From: "Kevin J. Sung" Date: Tue, 29 Oct 2024 21:37:40 -0400 Subject: [PATCH] remove deprecated classes from variational module --- python/ffsim/__init__.py | 6 - python/ffsim/qiskit/__init__.py | 6 - python/ffsim/qiskit/gates/__init__.py | 6 - python/ffsim/qiskit/gates/givens_ansatz.py | 92 +-- python/ffsim/qiskit/gates/ucj_operator.py | 84 -- python/ffsim/qiskit/sim.py | 24 +- python/ffsim/random/__init__.py | 2 - python/ffsim/random/random.py | 45 -- python/ffsim/variational/__init__.py | 6 +- python/ffsim/variational/givens.py | 70 -- python/ffsim/variational/ucj.py | 723 ------------------ .../python/qiskit/gates/givens_ansatz_test.py | 61 -- .../python/qiskit/gates/ucj_operator_test.py | 47 -- tests/python/qiskit/sampler_test.py | 30 +- tests/python/variational/givens_test.py | 54 -- tests/python/variational/ucj_test.py | 484 ------------ 16 files changed, 20 insertions(+), 1720 deletions(-) delete mode 100644 python/ffsim/qiskit/gates/ucj_operator.py delete mode 100644 python/ffsim/variational/ucj.py delete mode 100644 tests/python/qiskit/gates/ucj_operator_test.py delete mode 100644 tests/python/variational/ucj_test.py diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py index 75fa2673f..c0f8db381 100644 --- a/python/ffsim/__init__.py +++ b/python/ffsim/__init__.py @@ -92,12 +92,9 @@ ) from ffsim.variational import ( GivensAnsatzOp, - GivensAnsatzOperator, HopGateAnsatzOperator, NumNumAnsatzOpSpinBalanced, - RealUCJOperator, UCCSDOpRestrictedReal, - UCJOperator, UCJOpSpinBalanced, UCJOpSpinless, UCJOpSpinUnbalanced, @@ -112,13 +109,11 @@ "FermionAction", "FermionOperator", "GivensAnsatzOp", - "GivensAnsatzOperator", "HopGateAnsatzOperator", "MolecularData", "MolecularHamiltonian", "NumNumAnsatzOpSpinBalanced", "ProductStateSum", - "RealUCJOperator", "SingleFactorizedHamiltonian", "Spin", "StateVector", @@ -129,7 +124,6 @@ "SupportsLinearOperator", "SupportsTrace", "UCCSDOpRestrictedReal", - "UCJOperator", "UCJOpSpinBalanced", "UCJOpSpinUnbalanced", "UCJOpSpinless", diff --git a/python/ffsim/qiskit/__init__.py b/python/ffsim/qiskit/__init__.py index d3e34a37c..be69b6bfc 100644 --- a/python/ffsim/qiskit/__init__.py +++ b/python/ffsim/qiskit/__init__.py @@ -15,8 +15,6 @@ from ffsim.qiskit.gates import ( DiagCoulombEvolutionJW, DiagCoulombEvolutionSpinlessJW, - GivensAnsatzOperatorJW, - GivensAnsatzOperatorSpinlessJW, GivensAnsatzOpJW, GivensAnsatzOpSpinlessJW, NumNumAnsatzOpSpinBalancedJW, @@ -30,7 +28,6 @@ PrepareSlaterDeterminantSpinlessJW, SimulateTrotterDiagCoulombSplitOpJW, SimulateTrotterDoubleFactorizedJW, - UCJOperatorJW, UCJOpSpinBalancedJW, UCJOpSpinlessJW, UCJOpSpinUnbalancedJW, @@ -57,8 +54,6 @@ "FfsimSampler", "GivensAnsatzOpJW", "GivensAnsatzOpSpinlessJW", - "GivensAnsatzOperatorJW", - "GivensAnsatzOperatorSpinlessJW", "MergeOrbitalRotations", "NumNumAnsatzOpSpinBalancedJW", "NumOpSumEvolutionJW", @@ -72,7 +67,6 @@ "PrepareSlaterDeterminantSpinlessJW", "SimulateTrotterDiagCoulombSplitOpJW", "SimulateTrotterDoubleFactorizedJW", - "UCJOperatorJW", "UCJOpSpinBalancedJW", "UCJOpSpinUnbalancedJW", "UCJOpSpinlessJW", diff --git a/python/ffsim/qiskit/gates/__init__.py b/python/ffsim/qiskit/gates/__init__.py index cee382d71..faaeeada0 100644 --- a/python/ffsim/qiskit/gates/__init__.py +++ b/python/ffsim/qiskit/gates/__init__.py @@ -19,8 +19,6 @@ SimulateTrotterDoubleFactorizedJW, ) from ffsim.qiskit.gates.givens_ansatz import ( - GivensAnsatzOperatorJW, - GivensAnsatzOperatorSpinlessJW, GivensAnsatzOpJW, GivensAnsatzOpSpinlessJW, ) @@ -44,15 +42,12 @@ UCJOpSpinlessJW, UCJOpSpinUnbalancedJW, ) -from ffsim.qiskit.gates.ucj_operator import UCJOperatorJW __all__ = [ "DiagCoulombEvolutionJW", "DiagCoulombEvolutionSpinlessJW", "GivensAnsatzOpJW", "GivensAnsatzOpSpinlessJW", - "GivensAnsatzOperatorJW", - "GivensAnsatzOperatorSpinlessJW", "NumNumAnsatzOpSpinBalancedJW", "NumOpSumEvolutionJW", "NumOpSumEvolutionSpinlessJW", @@ -64,7 +59,6 @@ "PrepareSlaterDeterminantSpinlessJW", "SimulateTrotterDiagCoulombSplitOpJW", "SimulateTrotterDoubleFactorizedJW", - "UCJOperatorJW", "UCJOpSpinBalancedJW", "UCJOpSpinUnbalancedJW", "UCJOpSpinlessJW", diff --git a/python/ffsim/qiskit/gates/givens_ansatz.py b/python/ffsim/qiskit/gates/givens_ansatz.py index 4acecf5a0..287e119e3 100644 --- a/python/ffsim/qiskit/gates/givens_ansatz.py +++ b/python/ffsim/qiskit/gates/givens_ansatz.py @@ -12,7 +12,6 @@ from __future__ import annotations -import itertools import math from collections.abc import Iterator, Sequence @@ -25,97 +24,8 @@ Qubit, ) from qiskit.circuit.library import PhaseGate, XXPlusYYGate -from typing_extensions import deprecated -from ffsim.variational import GivensAnsatzOp, GivensAnsatzOperator - - -@deprecated("GivensAnsatzOperatorJW is deprecated. Use GivensAnsatzOpJW instead.") -class GivensAnsatzOperatorJW(Gate): - """Givens rotation ansatz operator under the Jordan-Wigner transformation. - - .. warning:: - This class is deprecated. Use :class:`ffsim.qiskit.GivensAnsatzOpJW` instead. - - See :class:`ffsim.GivensAnsatzOperator` for a description of this gate's unitary. - """ - - def __init__( - self, givens_ansatz_operator: GivensAnsatzOperator, *, label: str | None = None - ): - """Create a new Givens ansatz operator gate. - - Args: - givens_ansatz_operator: The Givens rotation ansatz operator. - label: The label of the gate. - """ - self.givens_ansatz_operator = givens_ansatz_operator - super().__init__( - "givens_ansatz_jw", 2 * givens_ansatz_operator.norb, [], label=label - ) - - def _define(self): - """Gate decomposition.""" - qubits = QuantumRegister(self.num_qubits) - circuit = QuantumCircuit(qubits, name=self.name) - norb = len(qubits) // 2 - alpha_qubits = qubits[:norb] - beta_qubits = qubits[norb:] - for instruction in _givens_ansatz_operator_jw( - alpha_qubits, self.givens_ansatz_operator - ): - circuit.append(instruction) - for instruction in _givens_ansatz_operator_jw( - beta_qubits, self.givens_ansatz_operator - ): - circuit.append(instruction) - self.definition = circuit - - -@deprecated( - "GivensAnsatzOperatorSpinlessJW is deprecated. " - "Use GivensAnsatzOpSpinlessJW instead." -) -class GivensAnsatzOperatorSpinlessJW(Gate): - """Givens rotation ansatz operator under the Jordan-Wigner transformation, spinless. - - Like :class:`GivensAnsatzOperatorJW` but only acts on a single spin species. - """ - - def __init__( - self, givens_ansatz_operator: GivensAnsatzOperator, *, label: str | None = None - ): - """Create a new Givens ansatz operator gate. - - Args: - givens_ansatz_operator: The Givens rotation ansatz operator. - label: The label of the gate. - """ - self.givens_ansatz_operator = givens_ansatz_operator - super().__init__( - "givens_ansatz_jw", givens_ansatz_operator.norb, [], label=label - ) - - def _define(self): - """Gate decomposition.""" - qubits = QuantumRegister(self.num_qubits) - self.definition = QuantumCircuit.from_instructions( - _givens_ansatz_operator_jw(qubits, self.givens_ansatz_operator), - qubits=qubits, - name=self.name, - ) - - -def _givens_ansatz_operator_jw( - qubits: Sequence[Qubit], givens_ansatz_operator: GivensAnsatzOperator -) -> Iterator[CircuitInstruction]: - for (i, j), theta in zip( - itertools.cycle(givens_ansatz_operator.interaction_pairs), - givens_ansatz_operator.thetas, - ): - yield CircuitInstruction( - XXPlusYYGate(2 * theta, -0.5 * math.pi), (qubits[i], qubits[j]) - ) +from ffsim.variational import GivensAnsatzOp class GivensAnsatzOpJW(Gate): diff --git a/python/ffsim/qiskit/gates/ucj_operator.py b/python/ffsim/qiskit/gates/ucj_operator.py deleted file mode 100644 index 03dfd0ff8..000000000 --- a/python/ffsim/qiskit/gates/ucj_operator.py +++ /dev/null @@ -1,84 +0,0 @@ -# (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. - -"""(Local) unitary cluster Jastrow ansatz gate.""" - -from __future__ import annotations - -from collections.abc import Iterator, Sequence - -from qiskit.circuit import ( - CircuitInstruction, - Gate, - QuantumCircuit, - QuantumRegister, - Qubit, -) -from typing_extensions import deprecated - -from ffsim.qiskit.gates.diag_coulomb import DiagCoulombEvolutionJW -from ffsim.qiskit.gates.orbital_rotation import OrbitalRotationJW -from ffsim.variational import UCJOperator - - -@deprecated("The UCJOperatorJW class is deprecated. Use UCJOpSpinBalancedJW instead.") -class UCJOperatorJW(Gate): - """Unitary cluster Jastrow operator under the Jordan-Wigner transformation. - - .. warning:: - This class is deprecated. Use :class:`ffsim.qiskit.UCJOpSpinBalancedJW` instead. - - See :class:`ffsim.UCJOperator` for a description of this gate's unitary. - - 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. - """ - - def __init__(self, ucj_operator: UCJOperator, *, label: str | None = None): - """Create a new unitary cluster Jastrow (UCJ) gate. - - Args: - ucj_operator: The UCJ operator. - label: The label of the gate. - """ - self.ucj_operator = ucj_operator - super().__init__("ucj_jw", 2 * ucj_operator.norb, [], label=label) - - def _define(self): - """Gate decomposition.""" - qubits = QuantumRegister(self.num_qubits) - self.definition = QuantumCircuit.from_instructions( - _ucj_jw(qubits, self.ucj_operator), qubits=qubits, name=self.name - ) - - -def _ucj_jw( - qubits: Sequence[Qubit], ucj_op: UCJOperator -) -> Iterator[CircuitInstruction]: - for mat, mat_alpha_beta, orbital_rotation in zip( - ucj_op.diag_coulomb_mats_alpha_alpha, - ucj_op.diag_coulomb_mats_alpha_beta, - ucj_op.orbital_rotations, - ): - yield CircuitInstruction( - OrbitalRotationJW(ucj_op.norb, orbital_rotation.T.conj()), qubits - ) - yield CircuitInstruction( - DiagCoulombEvolutionJW(ucj_op.norb, (mat, mat_alpha_beta, mat), -1.0), - qubits, - ) - yield CircuitInstruction( - OrbitalRotationJW(ucj_op.norb, orbital_rotation), qubits - ) - if ucj_op.final_orbital_rotation is not None: - yield CircuitInstruction( - OrbitalRotationJW(ucj_op.norb, ucj_op.final_orbital_rotation), qubits - ) diff --git a/python/ffsim/qiskit/sim.py b/python/ffsim/qiskit/sim.py index 7fd8dd316..2d0b521b3 100644 --- a/python/ffsim/qiskit/sim.py +++ b/python/ffsim/qiskit/sim.py @@ -21,8 +21,8 @@ from ffsim.qiskit.gates import ( DiagCoulombEvolutionJW, DiagCoulombEvolutionSpinlessJW, - GivensAnsatzOperatorJW, - GivensAnsatzOperatorSpinlessJW, + GivensAnsatzOpJW, + GivensAnsatzOpSpinlessJW, OrbitalRotationJW, OrbitalRotationSpinlessJW, PrepareHartreeFockJW, @@ -30,7 +30,6 @@ PrepareSlaterDeterminantJW, PrepareSlaterDeterminantSpinlessJW, SimulateTrotterDoubleFactorizedJW, - UCJOperatorJW, UCJOpSpinBalancedJW, UCJOpSpinlessJW, UCJOpSpinUnbalancedJW, @@ -136,14 +135,14 @@ def _evolve_state_vector_spinless( ) return states.StateVector(vec=vec, norb=norb, nelec=nelec) - if isinstance(op, GivensAnsatzOperatorSpinlessJW): + if isinstance(op, GivensAnsatzOpSpinlessJW): if not consecutive_sorted: raise ValueError( f"Gate of type '{op.__class__.__name__}' must be applied to " "consecutive qubits, in ascending order." ) vec = protocols.apply_unitary( - vec, op.givens_ansatz_operator, norb=norb, nelec=nelec, copy=False + vec, op.givens_ansatz_op, norb=norb, nelec=nelec, copy=False ) return states.StateVector(vec=vec, norb=norb, nelec=nelec) @@ -206,14 +205,14 @@ def _evolve_state_vector_spinful( ) return states.StateVector(vec=vec, norb=norb, nelec=nelec) - if isinstance(op, GivensAnsatzOperatorJW): + if isinstance(op, GivensAnsatzOpJW): if not consecutive_sorted: raise ValueError( f"Gate of type '{op.__class__.__name__}' must be applied to " "consecutive qubits, in ascending order." ) vec = protocols.apply_unitary( - vec, op.givens_ansatz_operator, norb=norb, nelec=nelec, copy=False + vec, op.givens_ansatz_op, norb=norb, nelec=nelec, copy=False ) return states.StateVector(vec=vec, norb=norb, nelec=nelec) @@ -263,17 +262,6 @@ def _evolve_state_vector_spinful( ) return states.StateVector(vec=vec, norb=norb, nelec=nelec) - if isinstance(op, UCJOperatorJW): - if not consecutive_sorted: - raise ValueError( - f"Gate of type '{op.__class__.__name__}' must be applied to " - "consecutive qubits, in ascending order." - ) - vec = protocols.apply_unitary( - vec, op.ucj_operator, norb=norb, nelec=nelec, copy=False - ) - return states.StateVector(vec=vec, norb=norb, nelec=nelec) - if isinstance(op, Barrier): return state_vector diff --git a/python/ffsim/random/__init__.py b/python/ffsim/random/__init__.py index 9ea9ed26c..4c0457d1f 100644 --- a/python/ffsim/random/__init__.py +++ b/python/ffsim/random/__init__.py @@ -29,7 +29,6 @@ random_ucj_op_spin_balanced, random_ucj_op_spin_unbalanced, random_ucj_op_spinless, - random_ucj_operator, random_unitary, ) @@ -49,7 +48,6 @@ "random_t2_amplitudes", "random_two_body_tensor", "random_uccsd_restricted", - "random_ucj_operator", "random_ucj_op_spin_balanced", "random_ucj_op_spin_unbalanced", "random_ucj_op_spinless", diff --git a/python/ffsim/random/random.py b/python/ffsim/random/random.py index 88697b275..d6ecf7055 100644 --- a/python/ffsim/random/random.py +++ b/python/ffsim/random/random.py @@ -300,51 +300,6 @@ def random_molecular_hamiltonian( ) -@deprecated( - "The random_ucj_operator function is deprecated. Use " - "random_ucj_operator_closed_shell or random_ucj_operator_open_shell instead." -) -def random_ucj_operator( - norb: int, - *, - n_reps: int = 1, - with_final_orbital_rotation: bool = False, - seed=None, -) -> variational.UCJOperator: - """Sample a random unitary cluster Jastrow (UCJ) operator. - - Args: - norb: The number of orbitals. - n_reps: The number of ansatz repetitions. - with_final_orbital_rotation: Whether to include a final orbital rotation - in the operator. - seed: A seed to initialize the pseudorandom number generator. - Should be a valid input to ``np.random.default_rng``. - - Returns: - The sampled UCJ operator. - """ - rng = np.random.default_rng(seed) - diag_coulomb_mats_alpha_alpha = np.stack( - [random_real_symmetric_matrix(norb, seed=rng) for _ in range(n_reps)] - ) - diag_coulomb_mats_alpha_beta = np.stack( - [random_real_symmetric_matrix(norb, seed=rng) for _ in range(n_reps)] - ) - orbital_rotations = np.stack( - [random_unitary(norb, seed=rng) for _ in range(n_reps)] - ) - final_orbital_rotation = None - if with_final_orbital_rotation: - final_orbital_rotation = random_unitary(norb, seed=rng) - return variational.UCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - final_orbital_rotation=final_orbital_rotation, - ) - - def random_uccsd_restricted( norb: int, nocc: int, diff --git a/python/ffsim/variational/__init__.py b/python/ffsim/variational/__init__.py index b330adcb9..33ca0324a 100644 --- a/python/ffsim/variational/__init__.py +++ b/python/ffsim/variational/__init__.py @@ -10,7 +10,7 @@ """Variational ansatzes.""" -from ffsim.variational.givens import GivensAnsatzOp, GivensAnsatzOperator +from ffsim.variational.givens import GivensAnsatzOp from ffsim.variational.hopgate import HopGateAnsatzOperator from ffsim.variational.multireference import ( multireference_state, @@ -18,19 +18,15 @@ ) from ffsim.variational.num_num import NumNumAnsatzOpSpinBalanced from ffsim.variational.uccsd import UCCSDOpRestrictedReal -from ffsim.variational.ucj import RealUCJOperator, UCJOperator from ffsim.variational.ucj_spin_balanced import UCJOpSpinBalanced from ffsim.variational.ucj_spin_unbalanced import UCJOpSpinUnbalanced from ffsim.variational.ucj_spinless import UCJOpSpinless __all__ = [ "GivensAnsatzOp", - "GivensAnsatzOperator", "HopGateAnsatzOperator", "NumNumAnsatzOpSpinBalanced", - "RealUCJOperator", "UCCSDOpRestrictedReal", - "UCJOperator", "UCJOpSpinBalanced", "UCJOpSpinUnbalanced", "UCJOpSpinless", diff --git a/python/ffsim/variational/givens.py b/python/ffsim/variational/givens.py index 43f0d6c3e..9ad5a4a33 100644 --- a/python/ffsim/variational/givens.py +++ b/python/ffsim/variational/givens.py @@ -19,82 +19,12 @@ from typing import cast import numpy as np -from scipy.linalg.blas import drot from scipy.linalg.lapack import zrot -from typing_extensions import deprecated from ffsim import linalg from ffsim.gates import apply_orbital_rotation -@dataclass(frozen=True) -@deprecated("GivensAnsatzOperator is deprecated. Use GivensAnsatzOp instead.") -class GivensAnsatzOperator: - """A Givens rotation ansatz operator. - - .. warning:: - This class is deprecated. Use :class:`ffsim.GivensAnsatzOp` instead. - - The Givens rotation ansatz consists of a sequence of `Givens rotations`_. - - Note that this ansatz does not implement any interactions between spin alpha and - spin beta orbitals. - - Attributes: - norb (int): The number of spatial orbitals. - interaction_pairs (list[tuple[int, int]]): The orbital pairs to apply the Givens - rotations to. - thetas (np.ndarray): The angles for the Givens rotations. - - .. _Givens rotations: ffsim.html#ffsim.apply_givens_rotation - """ - - norb: int - interaction_pairs: list[tuple[int, int]] - thetas: np.ndarray - - def _apply_unitary_( - self, vec: np.ndarray, norb: int, nelec: int | tuple[int, int], copy: bool - ) -> np.ndarray: - """Apply the operator to a vector.""" - return apply_orbital_rotation( - vec, self.to_orbital_rotation(), norb=norb, nelec=nelec, copy=copy - ) - - def to_parameters(self) -> np.ndarray: - """Convert the operator to a real-valued parameter vector.""" - num_params = len(self.thetas) - params = np.zeros(num_params) - params[: len(self.thetas)] = self.thetas - return params - - @staticmethod - def from_parameters( - params: np.ndarray, norb: int, interaction_pairs: list[tuple[int, int]] - ) -> GivensAnsatzOperator: - """Initialize the operator from a real-valued parameter vector. - - Args: - params: The real-valued parameter vector. - norb: The number of spatial orbitals. - interaction_pairs: The orbital pairs to apply the Givens rotation gates to. - """ - return GivensAnsatzOperator( - norb=norb, interaction_pairs=interaction_pairs, thetas=params - ) - - def to_orbital_rotation(self) -> np.ndarray: - orbital_rotation = np.eye(self.norb) - for (i, j), theta in zip(self.interaction_pairs[::-1], self.thetas[::-1]): - orbital_rotation[:, j], orbital_rotation[:, i] = drot( - orbital_rotation[:, j], - orbital_rotation[:, i], - math.cos(theta), - math.sin(theta), - ) - return orbital_rotation - - @dataclass(frozen=True) class GivensAnsatzOp: """A Givens rotation ansatz operator. diff --git a/python/ffsim/variational/ucj.py b/python/ffsim/variational/ucj.py deleted file mode 100644 index 1dfc70ab8..000000000 --- a/python/ffsim/variational/ucj.py +++ /dev/null @@ -1,723 +0,0 @@ -# (C) Copyright IBM 2023. -# -# 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. - -"""(Local) unitary cluster Jastrow ansatz.""" - -from __future__ import annotations - -import itertools -from dataclasses import dataclass -from typing import cast - -import numpy as np -import scipy.linalg -from opt_einsum import contract -from typing_extensions import deprecated - -from ffsim.gates import apply_diag_coulomb_evolution, apply_orbital_rotation -from ffsim.linalg import double_factorized_t2 -from ffsim.variational.util import ( - orbital_rotation_from_parameters, - orbital_rotation_to_parameters, -) - - -def _validate_diag_coulomb_indices(indices: list[tuple[int, int]] | None): - if indices is not None: - for i, j in indices: - if i > j: - raise ValueError( - "When specifying diagonal Coulomb indices, you must give only " - "upper trianglular indices. " - f"Got {(i, j)}, which is a lower triangular index." - ) - - -@dataclass(frozen=True) -@deprecated("The UCJOperator class is deprecated. Use UCJOpSpinBalanced instead.") -class UCJOperator: - r"""A unitary cluster Jastrow operator. - - .. warning:: - The UCJOperator class is deprecated. Use :class:`ffsim.UCJOpSpinBalanced` - instead. - - A unitary cluster Jastrow (UCJ) operator has the form - - .. math:: - - \prod_{k = 1}^L \mathcal{W}_k e^{i \mathcal{J}_k} \mathcal{W}_k^\dagger - - where each :math:`\mathcal{W_k}` is an orbital rotation and each :math:`\mathcal{J}` - is a diagonal Coulomb operator of the form - - .. math:: - - \mathcal{J} = \frac12\sum_{\sigma \tau, ij} - \mathbf{J}^{\sigma \tau}_{ij} n_{\sigma, i} n_{\tau, j}. - - In order that the operator commutes with the total spin Z operator, we enforce that - :math:`\mathbf{J}^{\alpha\alpha} = \mathbf{J}^{\beta\beta}` and - :math:`\mathbf{J}^{\alpha\beta} = \mathbf{J}^{\beta\alpha}`. As a result, we have - two sets of matrices for describing the diagonal Coulomb operators: - "alpha-alpha" matrices containing coefficients for terms involving the same spin, - and "alpha-beta" matrices containing coefficients for terms involving different - spins. - - To support variational optimization of the orbital basis, an optional final - orbital rotation can be included in the operator, to be performed at the end. - - Attributes: - diag_coulomb_mats_alpha_alpha (np.ndarray): The "alpha-alpha" diagonal Coulomb - matrices. - diag_coulomb_mats_alpha_beta (np.ndarray): The "alpha-beta" diagonal Coulomb - matrices. - orbital_rotations (np.ndarray): The orbital rotations. - final_orbital_rotation (np.ndarray): The optional final orbital rotation. - """ - - diag_coulomb_mats_alpha_alpha: np.ndarray - diag_coulomb_mats_alpha_beta: np.ndarray - orbital_rotations: np.ndarray - final_orbital_rotation: np.ndarray | None = None - - @property - def norb(self): - """The number of spatial orbitals.""" - return self.diag_coulomb_mats_alpha_alpha.shape[1] - - @property - def n_reps(self): - """The number of ansatz repetitions.""" - return self.diag_coulomb_mats_alpha_alpha.shape[0] - - @staticmethod - def n_params( - norb: int, - n_reps: int, - *, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - with_final_orbital_rotation: bool = False, - ) -> int: - """Return the number of parameters of an ansatz with given settings.""" - _validate_diag_coulomb_indices(alpha_alpha_indices) - _validate_diag_coulomb_indices(alpha_beta_indices) - n_params_aa = ( - norb * (norb + 1) // 2 - if alpha_alpha_indices is None - else len(alpha_alpha_indices) - ) - n_params_ab = ( - norb * (norb + 1) // 2 - if alpha_beta_indices is None - else len(alpha_beta_indices) - ) - return ( - n_reps * (norb**2 + n_params_aa + n_params_ab) - + with_final_orbital_rotation * norb**2 - ) - - @staticmethod - @deprecated("The UCJOperator class is deprecated. Use UCJOpSpinBalanced instead.") - def from_parameters( - params: np.ndarray, - *, - norb: int, - n_reps: int, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - with_final_orbital_rotation: bool = False, - ) -> UCJOperator: - r"""Initialize the UCJ operator from a real-valued parameter vector. - - Args: - params: The real-valued parameter vector. - norb: The number of spatial orbitals. - n_reps: The number of ansatz repetitions (:math:`L` from the docstring of - this class). - alpha_alpha_indices: Allowed indices for nonzero values of the "alpha-alpha" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - alpha_beta_indices: Allowed indices for nonzero values of the "alpha-beta" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - with_final_orbital_rotation: Whether to include a final orbital rotation - in the operator. - - Returns: - The UCJ operator constructed from the given parameters. - - Raises: - ValueError: alpha_alpha_indices contains lower triangular indices. - ValueError: alpha_beta_indices contains lower triangular indices. - """ - _validate_diag_coulomb_indices(alpha_alpha_indices) - _validate_diag_coulomb_indices(alpha_beta_indices) - return UCJOperator( - *_ucj_from_parameters( - params, - norb=norb, - n_reps=n_reps, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - with_final_orbital_rotation=with_final_orbital_rotation, - ) - ) - - def to_parameters( - self, - *, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - ) -> np.ndarray: - r"""Convert the UCJ operator to a real-valued parameter vector. - - If `alpha_alpha_indices` or `alpha_beta_indices` is specified, the returned - parameter vector will incorporate only the diagonal Coulomb matrix entries - corresponding to the given indices, so the original operator will not be - recoverable from the parameter vector. - - Args: - alpha_alpha_indices: Allowed indices for nonzero values of the "alpha-alpha" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - alpha_beta_indices: Allowed indices for nonzero values of the "alpha-beta" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - - Returns: - The real-valued parameter vector. - - Raises: - ValueError: alpha_alpha_indices contains lower triangular indices. - ValueError: alpha_beta_indices contains lower triangular indices. - """ - _validate_diag_coulomb_indices(alpha_alpha_indices) - _validate_diag_coulomb_indices(alpha_beta_indices) - return _ucj_to_parameters( - diag_coulomb_mats_alpha_alpha=self.diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=self.diag_coulomb_mats_alpha_beta, - orbital_rotations=self.orbital_rotations, - final_orbital_rotation=self.final_orbital_rotation, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - ) - - @staticmethod - @deprecated("The UCJOperator class is deprecated. Use UCJOpSpinBalanced instead.") - def from_t_amplitudes( - t2: np.ndarray, - *, - t1: np.ndarray | None = None, - n_reps: int | None = None, - tol: float = 1e-8, - ) -> UCJOperator: - """Initialize the UCJ operator from t2 (and optionally t1) amplitudes.""" - nocc, _, nvrt, _ = t2.shape - norb = nocc + nvrt - - diag_coulomb_mats, orbital_rotations = double_factorized_t2(t2, tol=tol) - diag_coulomb_mats = diag_coulomb_mats.reshape(-1, norb, norb)[:n_reps] - orbital_rotations = orbital_rotations.reshape(-1, norb, norb)[:n_reps] - - final_orbital_rotation = None - if t1 is not None: - final_orbital_rotation_generator = np.zeros((norb, norb), dtype=complex) - final_orbital_rotation_generator[:nocc, nocc:] = t1 - final_orbital_rotation_generator[nocc:, :nocc] = -t1.T - final_orbital_rotation = scipy.linalg.expm(final_orbital_rotation_generator) - - return UCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats, - orbital_rotations=orbital_rotations, - final_orbital_rotation=final_orbital_rotation, - ) - - def to_t_amplitudes( - self, nocc: int | None = None - ) -> np.ndarray | tuple[np.ndarray, np.ndarray]: - """Convert the UCJ operator to t2 (and possibly t1) amplitudes.""" - t2 = ( - 1j - * contract( - "kpq,kap,kip,kbq,kjq->ijab", - self.diag_coulomb_mats_alpha_alpha, - self.orbital_rotations, - self.orbital_rotations.conj(), - self.orbital_rotations, - self.orbital_rotations.conj(), - optimize="greedy", - )[:nocc, :nocc, nocc:, nocc:] - ) - - if self.final_orbital_rotation is None: - return t2 - - t1 = scipy.linalg.logm(self.final_orbital_rotation)[:nocc, nocc:] - return t2, t1 - - def _apply_unitary_( - self, vec: np.ndarray, norb: int, nelec: int | tuple[int, int], copy: bool - ) -> np.ndarray: - """Apply the operator to a vector.""" - if isinstance(nelec, int): - return NotImplemented - if copy: - vec = vec.copy() - for mat, mat_alpha_beta, orbital_rotation in zip( - self.diag_coulomb_mats_alpha_alpha, - self.diag_coulomb_mats_alpha_beta, - self.orbital_rotations, - ): - vec = apply_diag_coulomb_evolution( - vec, - mat=(mat, mat_alpha_beta, mat), - time=-1.0, - norb=norb, - nelec=nelec, - orbital_rotation=orbital_rotation, - copy=False, - ) - if self.final_orbital_rotation is not None: - vec = apply_orbital_rotation( - vec, - mat=self.final_orbital_rotation, - norb=norb, - nelec=nelec, - copy=False, - ) - return vec - - -@dataclass -@deprecated("The RealUCJOperator class is deprecated. Use UCJOpSpinBalanced instead.") -class RealUCJOperator: - r"""Real-valued unitary cluster Jastrow operator. - - .. warning:: - The RealUCJOperator class is deprecated. Use :class:`ffsim.UCJOpSpinBalanced` - instead. - - A real-valued unitary cluster Jastrow (UCJ) operator has the form - - .. math:: - - \prod_{k = 1}^L - \mathcal{W_k^*} e^{i \mathcal{-J}_k} \mathcal{W_k}^T - \mathcal{W_k} e^{i \mathcal{J}_k} \mathcal{W_k^\dagger} - - where each :math:`\mathcal{W_k}` is an orbital rotation and each :math:`\mathcal{J}` - is a diagonal Coulomb operator of the form - - .. math:: - - \mathcal{J} = \frac12\sum_{ij,\sigma \tau} - \mathbf{J}^{\sigma \tau}_{ij} n_{i,\sigma} n_{j,\tau}. - - In order that the operator commutes with the total spin Z operator, we enforce that - :math:`\mathbf{J}^{\alpha\alpha} = \mathbf{J}^{\beta\beta}` and - :math:`\mathbf{J}^{\alpha\beta} = \mathbf{J}^{\beta\alpha}`. As a result, we have - two sets of matrices for describing the diagonal Coulomb operators: - "alpha-alpha" matrices containing coefficients for terms involving the same spin, - and "alpha-beta" matrices containing coefficients for terms involving different - spins. - - To support variational optimization of the orbital basis, an optional final - orbital rotation can be included in the operator, to be performed at the end. - - Attributes: - diag_coulomb_mats_alpha_alpha: The "alpha-alpha" diagonal Coulomb matrices. - diag_coulomb_mats_alpha_beta: The "alpha-beta" diagonal Coulomb matrices. - orbital_rotations: The orbital rotations. - final_orbital_rotation: The optional final orbital rotation. - """ - - diag_coulomb_mats_alpha_alpha: np.ndarray - diag_coulomb_mats_alpha_beta: np.ndarray - orbital_rotations: np.ndarray - final_orbital_rotation: np.ndarray | None = None - - @property - def norb(self): - """The number of spatial orbitals.""" - return self.diag_coulomb_mats_alpha_alpha.shape[1] - - @property - def n_reps(self): - """The number of ansatz repetitions.""" - return self.diag_coulomb_mats_alpha_alpha.shape[0] - - @staticmethod - def n_params( - norb: int, - n_reps: int, - *, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - with_final_orbital_rotation: bool = False, - ) -> int: - """Return the number of parameters of an ansatz with given settings.""" - n_params_aa = ( - norb * (norb + 1) // 2 - if alpha_alpha_indices is None - else len(alpha_alpha_indices) - ) - n_params_ab = ( - norb * (norb + 1) // 2 - if alpha_beta_indices is None - else len(alpha_beta_indices) - ) - return ( - n_reps * (norb**2 + n_params_aa + n_params_ab) - + with_final_orbital_rotation * norb**2 - ) - - @staticmethod - @deprecated( - "The RealUCJOperator class is deprecated. Use UCJOpSpinBalanced instead." - ) - def from_parameters( - params: np.ndarray, - *, - norb: int, - n_reps: int, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - with_final_orbital_rotation: bool = False, - ) -> "RealUCJOperator": - r"""Initialize the real UCJ operator from a real-valued parameter vector. - - Args: - params: The real-valued parameter vector. - norb: The number of spatial orbitals. - n_reps: The number of ansatz repetitions (:math:`L` from the docstring of - this class). - alpha_alpha_indices: Allowed indices for nonzero values of the "alpha-alpha" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - alpha_beta_indices: Allowed indices for nonzero values of the "alpha-beta" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - with_final_orbital_rotation: Whether to include a final orbital rotation - in the operator. - - Returns: - The real UCJ operator constructed from the given parameters. - - Raises: - ValueError: alpha_alpha_indices contains lower triangular indices. - ValueError: alpha_beta_indices contains lower triangular indices. - """ - _validate_diag_coulomb_indices(alpha_alpha_indices) - _validate_diag_coulomb_indices(alpha_beta_indices) - return RealUCJOperator( - *_ucj_from_parameters( - params, - norb=norb, - n_reps=n_reps, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - with_final_orbital_rotation=with_final_orbital_rotation, - ) - ) - - def to_parameters( - self, - *, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - ) -> np.ndarray: - r"""Convert the real UCJ operator to a real-valued parameter vector. - - If `alpha_alpha_indices` or `alpha_beta_indices` is specified, the returned - parameter vector will incorporate only the diagonal Coulomb matrix entries - corresponding to the given indices, so the original operator will not be - recoverable from the parameter vector. - - Args: - alpha_alpha_indices: Allowed indices for nonzero values of the "alpha-alpha" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - alpha_beta_indices: Allowed indices for nonzero values of the "alpha-beta" - diagonal Coulomb matrices (see the docstring of this class). - If not specified, all matrix entries are allowed to be nonzero. - This list should contain only upper trianglular indices, i.e., - pairs :math:`(i, j)` where :math:`i \leq j`. Passing a list with - lower triangular indices will raise an error. - - Returns: - The real-valued parameter vector. - - Raises: - ValueError: alpha_alpha_indices contains lower triangular indices. - ValueError: alpha_beta_indices contains lower triangular indices. - """ - _validate_diag_coulomb_indices(alpha_alpha_indices) - _validate_diag_coulomb_indices(alpha_beta_indices) - return _ucj_to_parameters( - diag_coulomb_mats_alpha_alpha=self.diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=self.diag_coulomb_mats_alpha_beta, - orbital_rotations=self.orbital_rotations, - final_orbital_rotation=self.final_orbital_rotation, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - ) - - @staticmethod - @deprecated( - "The RealUCJOperator class is deprecated. Use UCJOpSpinBalanced instead." - ) - def from_t_amplitudes( - t2: np.ndarray, - *, - t1: np.ndarray | None = None, - n_reps: int | None = None, - tol: float = 1e-8, - ) -> "RealUCJOperator": - """Initialize the real UCJ operator from t2 (and optionally t1) amplitudes.""" - nocc, _, nvrt, _ = t2.shape - norb = nocc + nvrt - - diag_coulomb_mats, orbital_rotations = double_factorized_t2(t2, tol=tol) - diag_coulomb_mats = diag_coulomb_mats[:n_reps, 0] - orbital_rotations = orbital_rotations[:n_reps, 0] - - final_orbital_rotation = None - if t1 is not None: - final_orbital_rotation_generator = np.zeros((norb, norb), dtype=complex) - final_orbital_rotation_generator[:nocc, nocc:] = t1 - final_orbital_rotation_generator[nocc:, :nocc] = -t1.T - final_orbital_rotation = scipy.linalg.expm(final_orbital_rotation_generator) - - return RealUCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats, - orbital_rotations=orbital_rotations, - final_orbital_rotation=final_orbital_rotation, - ) - - def to_t_amplitudes( - self, nocc: int | None = None - ) -> np.ndarray | tuple[np.ndarray, np.ndarray]: - """Convert the UCJ operator to t2 (and possibly t1) amplitudes.""" - t2 = 1j * ( - contract( - "kpq,kap,kip,kbq,kjq->ijab", - self.diag_coulomb_mats_alpha_alpha, - self.orbital_rotations, - self.orbital_rotations.conj(), - self.orbital_rotations, - self.orbital_rotations.conj(), - optimize="greedy", - )[:nocc, :nocc, nocc:, nocc:] - + contract( - "kpq,kap,kip,kbq,kjq->ijab", - -self.diag_coulomb_mats_alpha_alpha, - self.orbital_rotations.conj(), - self.orbital_rotations, - self.orbital_rotations.conj(), - self.orbital_rotations, - optimize="greedy", - )[:nocc, :nocc, nocc:, nocc:] - ) - - if self.final_orbital_rotation is None: - return t2 - - t1 = scipy.linalg.logm(self.final_orbital_rotation)[:nocc, nocc:] - return t2, t1 - - def _apply_unitary_( - self, vec: np.ndarray, norb: int, nelec: int | tuple[int, int], copy: bool - ) -> np.ndarray: - """Apply the operator to a vector.""" - if isinstance(nelec, int): - return NotImplemented - if copy: - vec = vec.copy() - for mat, mat_alpha_beta, orbital_rotation in zip( - self.diag_coulomb_mats_alpha_alpha, - self.diag_coulomb_mats_alpha_beta, - self.orbital_rotations, - ): - vec = apply_diag_coulomb_evolution( - vec, - mat=(mat, mat_alpha_beta, mat), - time=-1.0, - norb=norb, - nelec=nelec, - orbital_rotation=orbital_rotation, - copy=False, - ) - vec = apply_diag_coulomb_evolution( - vec, - mat=(-mat, -mat_alpha_beta, -mat), - time=-1.0, - norb=norb, - nelec=nelec, - orbital_rotation=orbital_rotation.conj(), - copy=False, - ) - if self.final_orbital_rotation is not None: - vec = apply_orbital_rotation( - vec, - mat=self.final_orbital_rotation, - norb=norb, - nelec=nelec, - copy=False, - ) - return vec - - -def _ucj_from_parameters( - params: np.ndarray, - *, - norb: int, - n_reps: int, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, - with_final_orbital_rotation: bool = False, -) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]: - triu_indices = cast( - list[tuple[int, int]], - list(itertools.combinations_with_replacement(range(norb), 2)), - ) - if alpha_alpha_indices is None: - alpha_alpha_indices = triu_indices - if alpha_beta_indices is None: - alpha_beta_indices = triu_indices - diag_coulomb_mats_alpha_alpha = np.zeros((n_reps, norb, norb)) - diag_coulomb_mats_alpha_beta = np.zeros((n_reps, norb, norb)) - orbital_rotations = np.zeros((n_reps, norb, norb), dtype=complex) - index = 0 - for ( - orbital_rotation, - diag_coulomb_mat_alpha_alpha, - diag_coulomb_mat_alpha_beta, - ) in zip( - orbital_rotations, - diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta, - ): - # orbital rotation - n_params = norb**2 - orbital_rotation[:] = orbital_rotation_from_parameters( - params[index : index + n_params], norb - ) - index += n_params - # diag coulomb matrix, alpha-alpha - if alpha_alpha_indices: - n_params = len(alpha_alpha_indices) - rows, cols = zip(*alpha_alpha_indices) - vals = params[index : index + n_params] - diag_coulomb_mat_alpha_alpha[rows, cols] = vals - diag_coulomb_mat_alpha_alpha[cols, rows] = vals - index += n_params - # diag coulomb matrix, alpha-beta - if alpha_beta_indices: - n_params = len(alpha_beta_indices) - rows, cols = zip(*alpha_beta_indices) - vals = params[index : index + n_params] - diag_coulomb_mat_alpha_beta[rows, cols] = vals - diag_coulomb_mat_alpha_beta[cols, rows] = vals - index += n_params - # final orbital rotation - final_orbital_rotation = None - if with_final_orbital_rotation: - final_orbital_rotation = orbital_rotation_from_parameters(params[index:], norb) - return ( - diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta, - orbital_rotations, - final_orbital_rotation, - ) - - -def _ucj_to_parameters( - diag_coulomb_mats_alpha_alpha: np.ndarray, - diag_coulomb_mats_alpha_beta: np.ndarray, - orbital_rotations: np.ndarray, - final_orbital_rotation: np.ndarray | None = None, - *, - alpha_alpha_indices: list[tuple[int, int]] | None = None, - alpha_beta_indices: list[tuple[int, int]] | None = None, -) -> np.ndarray: - n_reps, norb, _ = diag_coulomb_mats_alpha_alpha.shape - triu_indices = cast( - list[tuple[int, int]], - list(itertools.combinations_with_replacement(range(norb), 2)), - ) - if alpha_alpha_indices is None: - alpha_alpha_indices = triu_indices - if alpha_beta_indices is None: - alpha_beta_indices = triu_indices - ntheta = n_reps * (len(alpha_alpha_indices) + len(alpha_beta_indices) + norb**2) - if final_orbital_rotation is not None: - ntheta += norb**2 - theta = np.zeros(ntheta) - index = 0 - for ( - orbital_rotation, - diag_coulomb_mat_alpha_alpha, - diag_coulomb_mat_alpha_beta, - ) in zip( - orbital_rotations, - diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta, - ): - # orbital rotation - n_params = norb**2 - theta[index : index + n_params] = orbital_rotation_to_parameters( - orbital_rotation - ) - index += n_params - # diag coulomb matrix, alpha-alpha - if alpha_alpha_indices: - n_params = len(alpha_alpha_indices) - theta[index : index + n_params] = diag_coulomb_mat_alpha_alpha[ - tuple(zip(*alpha_alpha_indices)) - ] - index += n_params - # diag coulomb matrix, alpha-beta - if alpha_beta_indices: - n_params = len(alpha_beta_indices) - theta[index : index + n_params] = diag_coulomb_mat_alpha_beta[ - tuple(zip(*alpha_beta_indices)) - ] - index += n_params - # final orbital rotation - if final_orbital_rotation is not None: - theta[index : index + norb**2] = orbital_rotation_to_parameters( - final_orbital_rotation - ) - return theta diff --git a/tests/python/qiskit/gates/givens_ansatz_test.py b/tests/python/qiskit/gates/givens_ansatz_test.py index 7ef4b91ef..477bc59a1 100644 --- a/tests/python/qiskit/gates/givens_ansatz_test.py +++ b/tests/python/qiskit/gates/givens_ansatz_test.py @@ -25,67 +25,6 @@ def brickwork(norb: int, n_layers: int): yield (j, j + 1) -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -@pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(5))) -def test_random_givens_ansatz_operator_spinful(norb: int, nelec: tuple[int, int]): - """Test random Givens rotation ansatz gives correct output state.""" - rng = np.random.default_rng() - dim = ffsim.dim(norb, nelec) - for _ in range(3): - n_layers = 2 * norb - interaction_pairs = list(brickwork(norb, n_layers)) - thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - givens_ansatz_op = ffsim.GivensAnsatzOperator(norb, interaction_pairs, thetas) - gate = ffsim.qiskit.GivensAnsatzOperatorJW(givens_ansatz_op) - - small_vec = ffsim.random.random_state_vector(dim, seed=rng) - big_vec = ffsim.qiskit.ffsim_vec_to_qiskit_vec( - small_vec, norb=norb, nelec=nelec - ) - - statevec = Statevector(big_vec).evolve(gate) - result = ffsim.qiskit.qiskit_vec_to_ffsim_vec( - np.array(statevec), norb=norb, nelec=nelec - ) - - expected = ffsim.apply_unitary( - small_vec, givens_ansatz_op, norb=norb, nelec=nelec - ) - - np.testing.assert_allclose(result, expected) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -@pytest.mark.parametrize("norb, nocc", ffsim.testing.generate_norb_nocc(range(5))) -def test_random_givens_ansatz_operator_spinless(norb: int, nocc: int): - """Test random spinless Givens rotation ansatz gives correct output state.""" - rng = np.random.default_rng() - nelec = (nocc, 0) - dim = ffsim.dim(norb, nelec) - for _ in range(3): - n_layers = 2 * norb - interaction_pairs = list(brickwork(norb, n_layers)) - thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - givens_ansatz_op = ffsim.GivensAnsatzOperator(norb, interaction_pairs, thetas) - gate = ffsim.qiskit.GivensAnsatzOperatorSpinlessJW(givens_ansatz_op) - - small_vec = ffsim.random.random_state_vector(dim, seed=rng) - big_vec = ffsim.qiskit.ffsim_vec_to_qiskit_vec( - small_vec, norb=norb, nelec=nelec - ) - - statevec = Statevector(big_vec).evolve(gate) - result = ffsim.qiskit.qiskit_vec_to_ffsim_vec( - np.array(statevec), norb=norb, nelec=nelec - ) - - expected = ffsim.apply_unitary( - small_vec, givens_ansatz_op, norb=norb, nelec=nelec - ) - - np.testing.assert_allclose(result, expected) - - @pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(5))) def test_random_spinful(norb: int, nelec: tuple[int, int]): """Test random Givens rotation ansatz gives correct output state.""" diff --git a/tests/python/qiskit/gates/ucj_operator_test.py b/tests/python/qiskit/gates/ucj_operator_test.py deleted file mode 100644 index 3f6d5e00a..000000000 --- a/tests/python/qiskit/gates/ucj_operator_test.py +++ /dev/null @@ -1,47 +0,0 @@ -# (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. - -"""Tests for unitary cluster Jastrow gate.""" - -from __future__ import annotations - -import numpy as np -import pytest -from qiskit.quantum_info import Statevector - -import ffsim - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -@pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(5))) -def test_random_ucj_operator(norb: int, nelec: tuple[int, int]): - """Test random UCJ gate gives correct output state.""" - rng = np.random.default_rng() - n_reps = 3 - dim = ffsim.dim(norb, nelec) - for _ in range(3): - ucj_op = ffsim.random.random_ucj_operator( - norb, n_reps=n_reps, with_final_orbital_rotation=True, seed=rng - ) - gate = ffsim.qiskit.UCJOperatorJW(ucj_op) - - small_vec = ffsim.random.random_state_vector(dim, seed=rng) - big_vec = ffsim.qiskit.ffsim_vec_to_qiskit_vec( - small_vec, norb=norb, nelec=nelec - ) - - statevec = Statevector(big_vec).evolve(gate) - result = ffsim.qiskit.qiskit_vec_to_ffsim_vec( - np.array(statevec), norb=norb, nelec=nelec - ) - - expected = ffsim.apply_unitary(small_vec, ucj_op, norb=norb, nelec=nelec) - - np.testing.assert_allclose(result, expected) diff --git a/tests/python/qiskit/sampler_test.py b/tests/python/qiskit/sampler_test.py index 180e08151..33fe17515 100644 --- a/tests/python/qiskit/sampler_test.py +++ b/tests/python/qiskit/sampler_test.py @@ -38,8 +38,6 @@ def _brickwork(norb: int, n_layers: int): yield (j, j + 1) -# TODO remove after removing UCJOperatorJW -@pytest.mark.filterwarnings("ignore::DeprecationWarning") @pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(1, 5))) def test_random_gates_spinful(norb: int, nelec: tuple[int, int]): """Test sampler with random gates.""" @@ -49,9 +47,6 @@ def test_random_gates_spinful(norb: int, nelec: tuple[int, int]): orbital_rotation = ffsim.random.random_unitary(norb, seed=rng) diag_coulomb_mat = ffsim.random.random_real_symmetric_matrix(norb, seed=rng) - ucj_op = ffsim.random.random_ucj_operator( - norb, n_reps=2, with_final_orbital_rotation=True, seed=rng - ) ucj_op_balanced = ffsim.random.random_ucj_op_spin_balanced( norb, n_reps=2, with_final_orbital_rotation=True, seed=rng ) @@ -66,7 +61,11 @@ def test_random_gates_spinful(norb: int, nelec: tuple[int, int]): ) interaction_pairs = list(_brickwork(norb, norb)) thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - givens_ansatz_op = ffsim.GivensAnsatzOperator(norb, interaction_pairs, thetas) + phis = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) + phase_angles = rng.uniform(-np.pi, np.pi, size=norb) + givens_ansatz_op = ffsim.GivensAnsatzOp( + norb, interaction_pairs, thetas, phis=phis, phase_angles=phase_angles + ) circuit = QuantumCircuit(qubits) circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits) @@ -74,8 +73,7 @@ def test_random_gates_spinful(norb: int, nelec: tuple[int, int]): circuit.append( ffsim.qiskit.DiagCoulombEvolutionJW(norb, diag_coulomb_mat, time=1.0), qubits ) - circuit.append(ffsim.qiskit.GivensAnsatzOperatorJW(givens_ansatz_op), qubits) - circuit.append(ffsim.qiskit.UCJOperatorJW(ucj_op), qubits) + circuit.append(ffsim.qiskit.GivensAnsatzOpJW(givens_ansatz_op), qubits) circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op_balanced), qubits) circuit.append(ffsim.qiskit.UCJOpSpinUnbalancedJW(ucj_op_unbalanced), qubits) circuit.append( @@ -111,8 +109,6 @@ def test_random_gates_spinful(norb: int, nelec: tuple[int, int]): assert np.sum(np.sqrt(exact_probs * empirical_probs)) > 0.999 -# TODO remove after removing UCJOperatorJW -@pytest.mark.filterwarnings("ignore::DeprecationWarning") @pytest.mark.parametrize("norb, nocc", ffsim.testing.generate_norb_nocc(range(1, 5))) def test_random_gates_spinless(norb: int, nocc: int): """Test sampler with random spinless gates.""" @@ -123,7 +119,11 @@ def test_random_gates_spinless(norb: int, nocc: int): orbital_rotation = ffsim.random.random_unitary(norb, seed=rng) interaction_pairs = list(_brickwork(norb, norb)) thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - givens_ansatz_op = ffsim.GivensAnsatzOperator(norb, interaction_pairs, thetas) + phis = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) + phase_angles = rng.uniform(-np.pi, np.pi, size=norb) + givens_ansatz_op = ffsim.GivensAnsatzOp( + norb, interaction_pairs, thetas, phis=phis, phase_angles=phase_angles + ) ucj_op = ffsim.random.random_ucj_op_spinless( norb, n_reps=2, with_final_orbital_rotation=True, seed=rng ) @@ -133,9 +133,7 @@ def test_random_gates_spinless(norb: int, nocc: int): circuit.append( ffsim.qiskit.OrbitalRotationSpinlessJW(norb, orbital_rotation), qubits ) - circuit.append( - ffsim.qiskit.GivensAnsatzOperatorSpinlessJW(givens_ansatz_op), qubits - ) + circuit.append(ffsim.qiskit.GivensAnsatzOpSpinlessJW(givens_ansatz_op), qubits) circuit.append(ffsim.qiskit.UCJOpSpinlessJW(ucj_op), qubits) circuit.measure_all() @@ -160,8 +158,6 @@ def test_random_gates_spinless(norb: int, nocc: int): assert np.sum(np.sqrt(exact_probs * empirical_probs)) > 0.999 -# TODO remove after removing UCJOperatorJW -@pytest.mark.filterwarnings("ignore::DeprecationWarning") @pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(1, 5))) def test_measure_subset_spinful(norb: int, nelec: tuple[int, int]): """Test measuring a subset of qubits.""" @@ -207,8 +203,6 @@ def test_measure_subset_spinful(norb: int, nelec: tuple[int, int]): assert _fidelity(ffsim_probs, qiskit_probs) > 0.99 -# TODO remove after removing UCJOperatorJW -@pytest.mark.filterwarnings("ignore::DeprecationWarning") @pytest.mark.parametrize("norb, nocc", ffsim.testing.generate_norb_nocc(range(1, 5))) def test_measure_subset_spinless(norb: int, nocc: int): """Test measuring a subset of qubits, spinless.""" diff --git a/tests/python/variational/givens_test.py b/tests/python/variational/givens_test.py index cfdd4efb1..fc50dd136 100644 --- a/tests/python/variational/givens_test.py +++ b/tests/python/variational/givens_test.py @@ -13,7 +13,6 @@ from __future__ import annotations import itertools -import math import numpy as np import pyscf @@ -24,59 +23,6 @@ import ffsim.linalg.givens -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_parameters_roundtrip(): - """Test converting to and back from parameters gives consistent results.""" - norb = 5 - n_reps = 2 - rng = np.random.default_rng() - - def ncycles(iterable, n): - "Returns the sequence elements n times" - return itertools.chain.from_iterable(itertools.repeat(tuple(iterable), n)) - - pairs = itertools.combinations(range(norb), 2) - interaction_pairs = list(ncycles(pairs, n_reps)) - thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - - operator = ffsim.GivensAnsatzOperator( - norb=norb, interaction_pairs=interaction_pairs, thetas=thetas - ) - assert len(operator.to_parameters()) == n_reps * math.comb(norb, 2) - roundtripped = ffsim.GivensAnsatzOperator.from_parameters( - operator.to_parameters(), norb=norb, interaction_pairs=interaction_pairs - ) - - np.testing.assert_allclose(roundtripped.thetas, operator.thetas) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -@pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nelec(range(5))) -def test_orbital_rotation(norb: int, nelec: tuple[int, int]): - n_reps = 2 - rng = np.random.default_rng() - - def ncycles(iterable, n): - "Returns the sequence elements n times" - return itertools.chain.from_iterable(itertools.repeat(tuple(iterable), n)) - - pairs = itertools.combinations(range(norb), 2) - interaction_pairs = list(ncycles(pairs, n_reps)) - thetas = rng.uniform(-np.pi, np.pi, size=len(interaction_pairs)) - - operator = ffsim.GivensAnsatzOperator( - norb=norb, interaction_pairs=interaction_pairs, thetas=thetas - ) - - vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng) - actual = ffsim.apply_orbital_rotation( - vec, operator.to_orbital_rotation(), norb=norb, nelec=nelec - ) - expected = ffsim.apply_unitary(vec, operator, norb=norb, nelec=nelec) - - np.testing.assert_allclose(actual, expected) - - def test_givens_parameters_roundtrip(): """Test converting to and back from parameters gives consistent results.""" norb = 5 diff --git a/tests/python/variational/ucj_test.py b/tests/python/variational/ucj_test.py deleted file mode 100644 index 2c3de272b..000000000 --- a/tests/python/variational/ucj_test.py +++ /dev/null @@ -1,484 +0,0 @@ -# (C) Copyright IBM 2023. -# -# 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. - -"""Tests for unitary cluster Jastrow ansatz.""" - -import itertools - -import numpy as np -import pyscf -import pytest -import scipy.linalg -from pyscf import cc - -import ffsim - - -def _exponentiate_t1(t1: np.ndarray, norb: int, nocc: int) -> np.ndarray: - generator = np.zeros((norb, norb), dtype=complex) - generator[:nocc, nocc:] = t1 - generator[nocc:, :nocc] = -t1.T - return scipy.linalg.expm(generator) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_n_params(): - for norb, n_reps, with_final_orbital_rotation in itertools.product( - [1, 2, 3], [1, 2, 3], [False, True] - ): - diag_coulomb_mats_alpha_alpha = np.zeros((n_reps, norb, norb)) - diag_coulomb_mats_alpha_beta = np.zeros((n_reps, norb, norb)) - orbital_rotations = np.stack([np.eye(norb) for _ in range(n_reps)]) - - final_orbital_rotation = np.eye(norb) - operator = ffsim.UCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - final_orbital_rotation=( - final_orbital_rotation if with_final_orbital_rotation else None - ), - ) - - actual = ffsim.UCJOperator.n_params( - norb, n_reps, with_final_orbital_rotation=with_final_orbital_rotation - ) - expected = len(operator.to_parameters()) - assert actual == expected - - alpha_alpha_indices = list( - itertools.combinations_with_replacement(range(norb), 2) - )[:norb] - alpha_beta_indices = list( - itertools.combinations_with_replacement(range(norb), 2) - )[:norb] - - actual = ffsim.UCJOperator.n_params( - norb, - n_reps, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - with_final_orbital_rotation=with_final_orbital_rotation, - ) - expected = len( - operator.to_parameters( - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - ) - ) - assert actual == expected - - with pytest.raises(ValueError, match="triangular"): - actual = ffsim.UCJOperator.n_params( - norb, - n_reps, - alpha_alpha_indices=[(1, 0)], - alpha_beta_indices=alpha_beta_indices, - ) - with pytest.raises(ValueError, match="triangular"): - actual = ffsim.UCJOperator.n_params( - norb, - n_reps, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=[(1, 0)], - ) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_parameters_roundtrip(): - norb = 5 - n_reps = 2 - diag_coulomb_mats_alpha_alpha = np.array( - [ffsim.random.random_real_symmetric_matrix(norb) for _ in range(n_reps)] - ) - diag_coulomb_mats_alpha_beta = np.array( - [ffsim.random.random_real_symmetric_matrix(norb) for _ in range(n_reps)] - ) - orbital_rotations = np.array( - [ffsim.random.random_unitary(norb) for _ in range(n_reps)] - ) - final_orbital_rotation = ffsim.random.random_unitary(norb) - operator = ffsim.UCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - final_orbital_rotation=final_orbital_rotation, - ) - roundtripped = ffsim.UCJOperator.from_parameters( - operator.to_parameters(), - norb=norb, - n_reps=n_reps, - with_final_orbital_rotation=True, - ) - np.testing.assert_allclose( - roundtripped.diag_coulomb_mats_alpha_alpha, - operator.diag_coulomb_mats_alpha_alpha, - ) - np.testing.assert_allclose( - roundtripped.diag_coulomb_mats_alpha_beta, - operator.diag_coulomb_mats_alpha_beta, - ) - np.testing.assert_allclose( - roundtripped.orbital_rotations, - operator.orbital_rotations, - ) - np.testing.assert_allclose( - roundtripped.final_orbital_rotation, - operator.final_orbital_rotation, - ) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_t_amplitudes_roundtrip(): - norb = 5 - nocc = 3 - - rng = np.random.default_rng() - - t2 = ffsim.random.random_t2_amplitudes(norb, nocc, dtype=float) - t1 = rng.standard_normal((nocc, norb - nocc)) - - operator = ffsim.UCJOperator.from_t_amplitudes(t2, t1=t1) - t2_roundtripped, t1_roundtripped = operator.to_t_amplitudes(nocc=nocc) - - np.testing.assert_allclose(t2_roundtripped, t2, atol=1e-12) - np.testing.assert_allclose( - _exponentiate_t1(t1_roundtripped, norb=norb, nocc=nocc), - _exponentiate_t1(t1, norb=norb, nocc=nocc), - atol=1e-12, - ) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_t_amplitudes(): - # Build an H2 molecule - mol = pyscf.gto.Mole() - mol.build( - atom=[["H", (0, 0, 0)], ["H", (0, 0, 1.8)]], - basis="sto-6g", - ) - hartree_fock = pyscf.scf.RHF(mol) - hartree_fock.kernel() - - # Get molecular data and molecular Hamiltonian (one- and two-body tensors) - mol_data = ffsim.MolecularData.from_scf(hartree_fock) - norb = mol_data.norb - nelec = mol_data.nelec - mol_hamiltonian = mol_data.hamiltonian - - # Get CCSD t2 amplitudes for initializing the ansatz - ccsd = cc.CCSD(hartree_fock) - _, _, t2 = ccsd.kernel() - - # Construct UCJ operator - n_reps = 2 - operator = ffsim.UCJOperator.from_t_amplitudes(t2, n_reps=n_reps) - - # Construct the Hartree-Fock state to use as the reference state - n_alpha, n_beta = nelec - reference_state = ffsim.slater_determinant( - norb=norb, occupied_orbitals=(range(n_alpha), range(n_beta)) - ) - - # Apply the operator to the reference state - ansatz_state = ffsim.apply_unitary( - reference_state, operator, norb=norb, nelec=nelec - ) - - # Compute the energy ⟨ψ|H|ψ⟩ of the ansatz state - hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) - energy = np.real(np.vdot(ansatz_state, hamiltonian @ ansatz_state)) - np.testing.assert_allclose(energy, -0.96962461) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_t_amplitudes_spin(): - """Test that initialization from CCSD amplitudes gives a singlet.""" - # Build an N2 molecule - mol = pyscf.gto.Mole() - mol.build( - atom=[["N", (0, 0, 0)], ["N", (0, 0, 1.5)]], - basis="sto-6g", - ) - hartree_fock = pyscf.scf.RHF(mol) - hartree_fock.kernel() - - # Get molecular data and molecular Hamiltonian (one- and two-body tensors) - mol_data = ffsim.MolecularData.from_scf(hartree_fock) - norb = mol_data.norb - nelec = mol_data.nelec - mol_hamiltonian = mol_data.hamiltonian - - # Get CCSD t2 amplitudes for initializing the ansatz - ccsd = cc.CCSD(hartree_fock) - _, _, t2 = ccsd.kernel() - - # Construct UCJ operator - n_reps = None - operator = ffsim.UCJOperator.from_t_amplitudes(t2, n_reps=n_reps) - - # Construct the Hartree-Fock state to use as the reference state - n_alpha, n_beta = nelec - reference_state = ffsim.slater_determinant( - norb=norb, occupied_orbitals=(range(n_alpha), range(n_beta)) - ) - spin_squared = ffsim.spin_square(reference_state, norb=norb, nelec=nelec) - np.testing.assert_allclose(spin_squared, 0) - - # Apply the operator to the reference state - ansatz_state = ffsim.apply_unitary( - reference_state, operator, norb=norb, nelec=nelec - ) - - # Compute the energy ⟨ψ|H|ψ⟩ of the ansatz state - hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) - energy = np.real(np.vdot(ansatz_state, hamiltonian @ ansatz_state)) - np.testing.assert_allclose(energy, -108.595692) - - # Compute the spin of the ansatz state - spin_squared = ffsim.spin_square(ansatz_state, norb=norb, nelec=nelec) - np.testing.assert_allclose(spin_squared, 0, atol=1e-12) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_n_params(): - for norb, n_reps, with_final_orbital_rotation in itertools.product( - [1, 2, 3], [1, 2, 3], [False, True] - ): - diag_coulomb_mats_alpha_alpha = np.zeros((n_reps, norb, norb)) - diag_coulomb_mats_alpha_beta = np.zeros((n_reps, norb, norb)) - orbital_rotations = np.stack([np.eye(norb) for _ in range(n_reps)]) - - final_orbital_rotation = np.eye(norb) - operator = ffsim.RealUCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - final_orbital_rotation=( - final_orbital_rotation if with_final_orbital_rotation else None - ), - ) - - actual = ffsim.RealUCJOperator.n_params( - norb, n_reps, with_final_orbital_rotation=with_final_orbital_rotation - ) - expected = len(operator.to_parameters()) - assert actual == expected - - alpha_alpha_indices = list( - itertools.combinations_with_replacement(range(norb), 2) - )[:norb] - alpha_beta_indices = list( - itertools.combinations_with_replacement(range(norb), 2) - )[:norb] - - actual = ffsim.RealUCJOperator.n_params( - norb, - n_reps, - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - with_final_orbital_rotation=with_final_orbital_rotation, - ) - expected = len( - operator.to_parameters( - alpha_alpha_indices=alpha_alpha_indices, - alpha_beta_indices=alpha_beta_indices, - ) - ) - assert actual == expected - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_parameters_roundtrip(): - norb = 5 - n_reps = 2 - diag_coulomb_mats_alpha_alpha = np.array( - [ffsim.random.random_real_symmetric_matrix(norb) for _ in range(n_reps)] - ) - diag_coulomb_mats_alpha_beta = np.array( - [ffsim.random.random_real_symmetric_matrix(norb) for _ in range(n_reps)] - ) - orbital_rotations = np.array( - [ffsim.random.random_unitary(norb) for _ in range(n_reps)] - ) - final_orbital_rotation = ffsim.random.random_unitary(norb) - operator = ffsim.RealUCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - final_orbital_rotation=final_orbital_rotation, - ) - roundtripped = ffsim.RealUCJOperator.from_parameters( - operator.to_parameters(), - norb=norb, - n_reps=n_reps, - with_final_orbital_rotation=True, - ) - np.testing.assert_allclose( - roundtripped.diag_coulomb_mats_alpha_alpha, - operator.diag_coulomb_mats_alpha_alpha, - ) - np.testing.assert_allclose( - roundtripped.diag_coulomb_mats_alpha_beta, - operator.diag_coulomb_mats_alpha_beta, - ) - np.testing.assert_allclose( - roundtripped.orbital_rotations, - operator.orbital_rotations, - ) - np.testing.assert_allclose( - roundtripped.final_orbital_rotation, - operator.final_orbital_rotation, - ) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_t_amplitudes_roundtrip(): - norb = 5 - nocc = 3 - - rng = np.random.default_rng() - - t2 = ffsim.random.random_t2_amplitudes(norb, nocc, dtype=float) - t1 = rng.standard_normal((nocc, norb - nocc)) - - operator = ffsim.RealUCJOperator.from_t_amplitudes(t2, t1=t1) - t2_roundtripped, t1_roundtripped = operator.to_t_amplitudes(nocc=nocc) - - np.testing.assert_allclose(t2_roundtripped, t2, atol=1e-12) - np.testing.assert_allclose( - _exponentiate_t1(t1_roundtripped, norb=norb, nocc=nocc), - _exponentiate_t1(t1, norb=norb, nocc=nocc), - atol=1e-12, - ) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_t_amplitudes(): - # Build an H2 molecule - mol = pyscf.gto.Mole() - mol.build( - atom=[["H", (0, 0, 0)], ["H", (0, 0, 1.8)]], - basis="sto-6g", - ) - hartree_fock = pyscf.scf.RHF(mol) - hartree_fock.kernel() - - # Get molecular data and molecular Hamiltonian (one- and two-body tensors) - mol_data = ffsim.MolecularData.from_scf(hartree_fock) - norb = mol_data.norb - nelec = mol_data.nelec - mol_hamiltonian = mol_data.hamiltonian - - # Get CCSD t2 amplitudes for initializing the ansatz - ccsd = cc.CCSD(hartree_fock) - _, _, t2 = ccsd.kernel() - - # Construct UCJ operator - n_reps = 1 - operator = ffsim.RealUCJOperator.from_t_amplitudes(t2, n_reps=n_reps) - - # Construct the Hartree-Fock state to use as the reference state - n_alpha, n_beta = nelec - reference_state = ffsim.slater_determinant( - norb=norb, occupied_orbitals=(range(n_alpha), range(n_beta)) - ) - - # Apply the operator to the reference state - ansatz_state = ffsim.apply_unitary( - reference_state, operator, norb=norb, nelec=nelec - ) - - # Compute the energy ⟨ψ|H|ψ⟩ of the ansatz state - hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) - energy = np.real(np.vdot(ansatz_state, hamiltonian @ ansatz_state)) - np.testing.assert_allclose(energy, -0.96962461) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_t_amplitudes_spin(): - """Test that initialization from CCSD amplitudes gives a singlet.""" - # Build an N2 molecule - mol = pyscf.gto.Mole() - mol.build( - atom=[["N", (0, 0, 0)], ["N", (0, 0, 1.5)]], - basis="sto-6g", - ) - hartree_fock = pyscf.scf.RHF(mol) - hartree_fock.kernel() - - # Get molecular data and molecular Hamiltonian (one- and two-body tensors) - mol_data = ffsim.MolecularData.from_scf(hartree_fock) - norb = mol_data.norb - nelec = mol_data.nelec - mol_hamiltonian = mol_data.hamiltonian - - # Get CCSD t2 amplitudes for initializing the ansatz - ccsd = cc.CCSD(hartree_fock) - _, _, t2 = ccsd.kernel() - - # Construct UCJ operator - n_reps = None - operator = ffsim.RealUCJOperator.from_t_amplitudes(t2, n_reps=n_reps) - - # Construct the Hartree-Fock state to use as the reference state - n_alpha, n_beta = nelec - reference_state = ffsim.slater_determinant( - norb=norb, occupied_orbitals=(range(n_alpha), range(n_beta)) - ) - spin_squared = ffsim.spin_square(reference_state, norb=norb, nelec=nelec) - np.testing.assert_allclose(spin_squared, 0) - - # Apply the operator to the reference state - ansatz_state = ffsim.apply_unitary( - reference_state, operator, norb=norb, nelec=nelec - ) - - # Compute the energy ⟨ψ|H|ψ⟩ of the ansatz state - hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) - energy = np.real(np.vdot(ansatz_state, hamiltonian @ ansatz_state)) - np.testing.assert_allclose(energy, -108.595692) - - # Compute the spin of the ansatz state - spin_squared = ffsim.spin_square(ansatz_state, norb=norb, nelec=nelec) - np.testing.assert_allclose(spin_squared, 0, atol=1e-12) - - -@pytest.mark.filterwarnings("ignore::DeprecationWarning") -def test_real_ucj_preserves_real(): - """Test that the real-valued UCJ ansatz preserves reality of t2 amplitudes.""" - norb = 5 - - rng = np.random.default_rng() - - n_reps = 2 - diag_coulomb_mats_alpha_alpha = np.stack( - [ - ffsim.random.random_real_symmetric_matrix(norb, seed=rng) - for _ in range(n_reps) - ] - ) - diag_coulomb_mats_alpha_beta = np.stack( - [ - ffsim.random.random_real_symmetric_matrix(norb, seed=rng) - for _ in range(n_reps) - ] - ) - orbital_rotations = np.stack( - [ffsim.random.random_unitary(norb, seed=rng) for _ in range(n_reps)] - ) - operator = ffsim.RealUCJOperator( - diag_coulomb_mats_alpha_alpha=diag_coulomb_mats_alpha_alpha, - diag_coulomb_mats_alpha_beta=diag_coulomb_mats_alpha_beta, - orbital_rotations=orbital_rotations, - ) - - t2 = operator.to_t_amplitudes() - np.testing.assert_allclose(np.imag(t2), 0)