diff --git a/pyproject.toml b/pyproject.toml index e1945e67c..96b12b79f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ dependencies = [ "numpy", "opt_einsum", "orjson", - "pyscf >= 2.4", + "pyscf >= 2.7", "qiskit >= 1.1", "scipy", "typing-extensions", diff --git a/python/ffsim/hamiltonians/molecular_hamiltonian.py b/python/ffsim/hamiltonians/molecular_hamiltonian.py index 806d6c5f3..0b6d8a721 100644 --- a/python/ffsim/hamiltonians/molecular_hamiltonian.py +++ b/python/ffsim/hamiltonians/molecular_hamiltonian.py @@ -18,7 +18,7 @@ import pyscf.ao2mo import pyscf.tools from opt_einsum import contract -from pyscf.fci.direct_nosym import absorb_h1e, contract_1e, contract_2e, make_hdiag +from pyscf.fci.direct_nosym import absorb_h1e, contract_2e, make_hdiag from scipy.sparse.linalg import LinearOperator from typing_extensions import deprecated @@ -122,34 +122,22 @@ def rotated(self, orbital_rotation: np.ndarray) -> MolecularHamiltonian: def _linear_operator_(self, norb: int, nelec: tuple[int, int]) -> LinearOperator: """Return a SciPy LinearOperator representing the object.""" - if np.iscomplexobj(self.two_body_tensor): - raise NotImplementedError( - "This Hamiltonian has a complex-valued two-body tensor. " - "LinearOperator support for complex two-body tensors is not yet " - "implemented. See https://github.com/qiskit-community/ffsim/issues/81." - ) - n_alpha, n_beta = nelec linkstr_index_a = gen_linkstr_index(range(norb), n_alpha) linkstr_index_b = gen_linkstr_index(range(norb), n_beta) link_index = (linkstr_index_a, linkstr_index_b) + two_body = absorb_h1e( + self.one_body_tensor, self.two_body_tensor, norb, nelec, 0.5 + ) - if np.iscomplexobj(self.one_body_tensor): - return _linear_operator_complex( - one_body_tensor=self.one_body_tensor, - two_body_tensor=self.two_body_tensor, - constant=self.constant, - norb=norb, - nelec=nelec, - link_index=link_index, - ) - return _linear_operator_real( - one_body_tensor=self.one_body_tensor, - two_body_tensor=self.two_body_tensor, - constant=self.constant, - norb=norb, - nelec=nelec, - link_index=link_index, + def matvec(vec: np.ndarray): + result = self.constant * vec.astype(complex, copy=False) + result += contract_2e(two_body, vec, norb, nelec, link_index=link_index) + return result + + dim_ = dim(norb, nelec) + return LinearOperator( + shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex ) def _diag_(self, norb: int, nelec: tuple[int, int]) -> np.ndarray: @@ -203,57 +191,3 @@ def _approx_eq_(self, other, rtol: float, atol: float) -> bool: return False return True return NotImplemented - - -def _linear_operator_complex( - one_body_tensor: np.ndarray, - two_body_tensor: np.ndarray, - constant: float, - norb: int, - nelec: tuple[int, int], - link_index: tuple[np.ndarray, np.ndarray], -): - two_body = absorb_h1e(one_body_tensor.real, two_body_tensor, norb, nelec, 0.5) - dim_ = dim(norb, nelec) - - def matvec(vec: np.ndarray): - result = constant * vec.astype(complex, copy=False) - result += 1j * contract_1e( - one_body_tensor.imag, vec.real, norb, nelec, link_index=link_index - ) - result -= contract_1e( - one_body_tensor.imag, vec.imag, norb, nelec, link_index=link_index - ) - result += contract_2e(two_body, vec.real, norb, nelec, link_index=link_index) - result += 1j * contract_2e( - two_body, vec.imag, norb, nelec, link_index=link_index - ) - return result - - return LinearOperator( - shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex - ) - - -def _linear_operator_real( - one_body_tensor: np.ndarray, - two_body_tensor: np.ndarray, - constant: float, - norb: int, - nelec: tuple[int, int], - link_index: tuple[np.ndarray, np.ndarray], -): - two_body = absorb_h1e(one_body_tensor, two_body_tensor, norb, nelec, 0.5) - dim_ = dim(norb, nelec) - - def matvec(vec: np.ndarray): - result = constant * vec.astype(complex, copy=False) - result += contract_2e(two_body, vec.real, norb, nelec, link_index=link_index) - result += 1j * contract_2e( - two_body, vec.imag, norb, nelec, link_index=link_index - ) - return result - - return LinearOperator( - shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex - ) diff --git a/tests/python/hamiltonians/molecular_hamiltonian_test.py b/tests/python/hamiltonians/molecular_hamiltonian_test.py index 1b56bdf2d..3b594970a 100644 --- a/tests/python/hamiltonians/molecular_hamiltonian_test.py +++ b/tests/python/hamiltonians/molecular_hamiltonian_test.py @@ -106,8 +106,7 @@ def test_fermion_operator(norb: int, nelec: tuple[int, int]): rng = np.random.default_rng() one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng) - # TODO remove dtype=float after adding support for complex - two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float) + two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng) constant = rng.standard_normal() mol_hamiltonian = ffsim.MolecularHamiltonian( one_body_tensor, two_body_tensor, constant=constant @@ -132,8 +131,7 @@ def test_rotated(): # generate a random molecular Hamiltonian one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng) - # TODO remove dtype=float after adding support for complex - two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float) + two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng) constant = rng.standard_normal() mol_hamiltonian = ffsim.MolecularHamiltonian( one_body_tensor, two_body_tensor, constant=constant diff --git a/tests/python/hamiltonians/single_factorized_hamiltonian_test.py b/tests/python/hamiltonians/single_factorized_hamiltonian_test.py index cbdae51e1..761a03007 100644 --- a/tests/python/hamiltonians/single_factorized_hamiltonian_test.py +++ b/tests/python/hamiltonians/single_factorized_hamiltonian_test.py @@ -33,7 +33,6 @@ def test_linear_operator(norb: int, nelec: tuple[int, int], cholesky: bool): """Test linear operator.""" rng = np.random.default_rng(2474) - # TODO remove dtype=float once complex two-body is supported mol_hamiltonian = ffsim.random.random_molecular_hamiltonian( norb, seed=rng, dtype=float ) @@ -58,7 +57,6 @@ def test_reduced_matrix_product_states(norb: int, nelec: tuple[int, int]): """Test computing reduced matrix on product states.""" rng = np.random.default_rng(7869) - # TODO remove dtype=float once complex two-body is supported mol_hamiltonian = ffsim.random.random_molecular_hamiltonian( norb, seed=rng, dtype=float ) @@ -104,7 +102,6 @@ def test_expectation_product_state_slater_determinant( """Test computing expectation value on Slater determinant product state.""" rng = np.random.default_rng(3400) - # TODO remove dtype=float once complex two-body is supported mol_hamiltonian = ffsim.random.random_molecular_hamiltonian( norb, seed=rng, dtype=float ) @@ -137,7 +134,6 @@ def test_expectation_product_state(norb: int, nelec: tuple[int, int]): """Test computing expectation value on product state.""" rng = np.random.default_rng(6775) - # TODO remove dtype=float once complex two-body is supported mol_hamiltonian = ffsim.random.random_molecular_hamiltonian( norb, seed=rng, dtype=float )