diff --git a/requirements.txt b/requirements.txt index 1044aaf56..1be5fc3fc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy>=1.11.0 +numpy>=1.13.0 scipy>=1.1.0 h5py>=2.8 future diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py index 8b3cd128a..b74da7b5e 100644 --- a/src/openfermion/ops/_quadratic_hamiltonian.py +++ b/src/openfermion/ops/_quadratic_hamiltonian.py @@ -12,6 +12,8 @@ """Hamiltonians that are quadratic in the fermionic ladder operators.""" +import warnings + import numpy from scipy.linalg import schur @@ -125,52 +127,12 @@ def add_chemical_potential(self, chemical_potential): numpy.eye(self.n_qubits)) self.chemical_potential += chemical_potential - def orbital_energies(self, non_negative=False): - r"""Return the orbital energies. - - Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian - of the form - - .. math:: - - \sum_{j} \varepsilon_j b^\dagger_j b_j + \text{constant}. - - We call the :math:`\varepsilon_j` the orbital energies. - The eigenvalues of the Hamiltonian are sums of subsets of the - orbital energies (up to the additive constant). - - Args: - non_negative(bool): If True, always return a list of orbital - energies that are non-negative. This option is ignored if - the Hamiltonian does not conserve particle number, in which - case the returned orbital energies are always non-negative. - - Returns - ------- - orbital_energies(ndarray) - A one-dimensional array containing the :math:`\varepsilon_j` - constant(float) - The constant - """ - if self.conserves_particle_number and not non_negative: - hermitian_matrix = self.combined_hermitian_part - orbital_energies, _ = numpy.linalg.eigh( - hermitian_matrix) - constant = self.constant - else: - majorana_matrix, majorana_constant = self.majorana_form() - canonical, _ = antisymmetric_canonical_form( - majorana_matrix) - orbital_energies = canonical[ - range(self.n_qubits), range(self.n_qubits, 2 * self.n_qubits)] - constant = -0.5 * numpy.sum(orbital_energies) + majorana_constant - - return orbital_energies, constant - def ground_energy(self): """Return the ground energy.""" - _, constant = self.orbital_energies(non_negative=True) - return constant + orbital_energies, _, constant = ( + self.diagonalizing_bogoliubov_transform()) + return numpy.sum(orbital_energies[ + numpy.where(orbital_energies < 0.0)[0]]) + constant def majorana_form(self): r"""Return the Majorana represention of the Hamiltonian. @@ -324,10 +286,30 @@ def index_map(i): spin_indices = [index_map(i) for i in range(n_sites)] matrix = self.combined_hermitian_part[ numpy.ix_(spin_indices, spin_indices)] + orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh( + matrix) else: matrix = self.combined_hermitian_part - orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh(matrix) + if _is_spin_block_diagonal(matrix): + up_block = matrix[:n_modes//2, :n_modes//2] + down_block = matrix[n_modes//2:, n_modes//2:] + + up_orbital_energies, up_diagonalizing_unitary_T = ( + numpy.linalg.eigh(up_block)) + down_orbital_energies, down_diagonalizing_unitary_T = ( + numpy.linalg.eigh(down_block)) + + orbital_energies = numpy.concatenate( + (up_orbital_energies, down_orbital_energies)) + diagonalizing_unitary_T = numpy.block([ + [up_diagonalizing_unitary_T, + numpy.zeros((n_modes//2, n_modes//2))], + [numpy.zeros((n_modes//2, n_modes//2)), + down_diagonalizing_unitary_T]]) + else: + orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh( + matrix) return orbital_energies, diagonalizing_unitary_T.T, self.constant @@ -417,6 +399,44 @@ def diagonalizing_circuit(self): return circuit_description + # DEPRECATED FUNCTIONS + # ==================== + def orbital_energies(self, non_negative=False): + r"""Return the orbital energies. + + Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian + of the form + + .. math:: + + \sum_{j} \varepsilon_j b^\dagger_j b_j + \text{constant}. + + We call the :math:`\varepsilon_j` the orbital energies. + The eigenvalues of the Hamiltonian are sums of subsets of the + orbital energies (up to the additive constant). + + Args: + non_negative(bool): If True, always return a list of orbital + energies that are non-negative. This option is ignored if + the Hamiltonian does not conserve particle number, in which + case the returned orbital energies are always non-negative. + + Returns + ------- + orbital_energies(ndarray) + A one-dimensional array containing the :math:`\varepsilon_j` + constant(float) + The constant + """ + warnings.warn('The method `orbital_energies` is deprecated. ' + 'Use the method `diagonalizing_bogoliubov_transform` ' + 'instead.', DeprecationWarning) + + orbital_energies, _, constant = ( + self.diagonalizing_bogoliubov_transform()) + + return orbital_energies, constant + def antisymmetric_canonical_form(antisymmetric_matrix): """Compute the canonical form of an antisymmetric matrix. @@ -495,3 +515,13 @@ def antisymmetric_canonical_form(antisymmetric_matrix): swap_rows(diagonal, i, arg_min) return canonical, orthogonal.T + + +def _is_spin_block_diagonal(matrix) -> bool: + n = matrix.shape[0] + if n % 2: + return False + max_upper_right = numpy.max(numpy.abs(matrix[:n//2, n//2:])) + max_lower_left = numpy.max(numpy.abs(matrix[n//2:, :n//2])) + return (numpy.isclose(max_upper_right, 0.0) + and numpy.isclose(max_lower_left, 0.0)) diff --git a/src/openfermion/ops/_quadratic_hamiltonian_test.py b/src/openfermion/ops/_quadratic_hamiltonian_test.py index d641967f4..385b47693 100644 --- a/src/openfermion/ops/_quadratic_hamiltonian_test.py +++ b/src/openfermion/ops/_quadratic_hamiltonian_test.py @@ -239,6 +239,29 @@ def test_diagonalizing_bogoliubov_transform_particle_conserving(self): quad_ham = get_quadratic_hamiltonian( reorder(get_fermion_operator(quad_ham), up_then_down)) + orbital_energies, transformation_matrix, _ = ( + quad_ham.diagonalizing_bogoliubov_transform() + ) + max_upper_right = numpy.max(numpy.abs(transformation_matrix[:5, 5:])) + max_lower_left = numpy.max(numpy.abs(transformation_matrix[5:, :5])) + + numpy.testing.assert_allclose( + orbital_energies[:5], orbital_energies[5:]) + numpy.testing.assert_allclose( + transformation_matrix.dot( + quad_ham.combined_hermitian_part.T.dot( + transformation_matrix.T.conj())), + numpy.diag(orbital_energies), + atol=1e-7) + numpy.testing.assert_allclose(max_upper_right, 0.0) + numpy.testing.assert_allclose(max_lower_left, 0.0) + + # Specific spin sector + quad_ham = random_quadratic_hamiltonian( + 5, conserves_particle_number=True, expand_spin=True) + quad_ham = get_quadratic_hamiltonian( + reorder(get_fermion_operator(quad_ham), up_then_down)) + for spin_sector in range(2): orbital_energies, transformation_matrix, _ = ( quad_ham.diagonalizing_bogoliubov_transform( diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py index 866344d86..31a9f862b 100644 --- a/src/openfermion/utils/_slater_determinants.py +++ b/src/openfermion/utils/_slater_determinants.py @@ -105,8 +105,7 @@ def gaussian_state_preparation_circuit( if occupied_orbitals is None: # The ground state is desired, so we fill the orbitals that have # negative energy - num_negative_energies = numpy.count_nonzero(orbital_energies < 0.0) - occupied_orbitals = range(num_negative_energies) + occupied_orbitals = numpy.where(orbital_energies < 0.0)[0] # Get the unitary rows which represent the Slater determinant slater_determinant_matrix = transformation_matrix[occupied_orbitals]