diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py
index c0f8db381..32181000a 100644
--- a/python/ffsim/__init__.py
+++ b/python/ffsim/__init__.py
@@ -71,19 +71,14 @@
     expectation_one_body_power,
     expectation_one_body_product,
     hartree_fock_state,
-    indices_to_strings,
-    one_hot,
-    rdm,
     rdms,
     sample_slater_determinant,
     sample_state_vector,
     slater_determinant,
     slater_determinant_amplitudes,
-    slater_determinant_rdm,
     slater_determinant_rdms,
     spin_square,
     strings_to_addresses,
-    strings_to_indices,
 )
 from ffsim.trotter import (
     simulate_qdrift_double_factorized,
@@ -158,18 +153,15 @@
     "fermi_hubbard_2d",
     "fermion_operator",
     "hartree_fock_state",
-    "indices_to_strings",
     "init_cache",
     "linalg",
     "linear_operator",
     "multireference_state",
     "multireference_state_prod",
     "number_operator",
-    "one_hot",
     "optimize",
     "qiskit",
     "random",
-    "rdm",
     "rdms",
     "sample_slater_determinant",
     "sample_state_vector",
@@ -178,11 +170,9 @@
     "simulate_trotter_double_factorized",
     "slater_determinant",
     "slater_determinant_amplitudes",
-    "slater_determinant_rdm",
     "slater_determinant_rdms",
     "spin_square",
     "strings_to_addresses",
-    "strings_to_indices",
     "testing",
     "trace",
 ]
diff --git a/python/ffsim/hamiltonians/molecular_hamiltonian.py b/python/ffsim/hamiltonians/molecular_hamiltonian.py
index 0b6d8a721..821223f69 100644
--- a/python/ffsim/hamiltonians/molecular_hamiltonian.py
+++ b/python/ffsim/hamiltonians/molecular_hamiltonian.py
@@ -12,15 +12,11 @@
 
 import dataclasses
 import itertools
-import os
 
 import numpy as np
-import pyscf.ao2mo
-import pyscf.tools
 from opt_einsum import contract
 from pyscf.fci.direct_nosym import absorb_h1e, contract_2e, make_hdiag
 from scipy.sparse.linalg import LinearOperator
-from typing_extensions import deprecated
 
 from ffsim.cistring import gen_linkstr_index
 from ffsim.operators import FermionOperator, cre_a, cre_b, des_a, des_b
@@ -53,29 +49,6 @@ class MolecularHamiltonian:
     two_body_tensor: np.ndarray
     constant: float = 0.0
 
-    @staticmethod
-    @deprecated(
-        "The MolecularHamiltonian.from_fcidump method is deprecated. "
-        "Instead, use MolecularData.from_fcidump and then access the `hamiltonian` "
-        "attribute of the returned MolecularData."
-    )
-    def from_fcidump(file: str | bytes | os.PathLike) -> MolecularHamiltonian:
-        """Initialize a MolecularHamiltonian from an FCIDUMP file.
-
-        .. warning::
-            This function is deprecated. Instead, use MolecularData.from_fcidump and
-            then access the `hamiltonian` attribute of the returned MolecularData.
-
-        Args:
-            file: The FCIDUMP file path.
-        """
-        data = pyscf.tools.fcidump.read(file, verbose=False)
-        return MolecularHamiltonian(
-            one_body_tensor=data["H1"],
-            two_body_tensor=pyscf.ao2mo.restore(1, data["H2"], data["NORB"]),
-            constant=data["ECORE"],
-        )
-
     @property
     def norb(self) -> int:
         """The number of spatial orbitals."""
diff --git a/python/ffsim/molecular_data.py b/python/ffsim/molecular_data.py
index 5a30aa051..23df81382 100644
--- a/python/ffsim/molecular_data.py
+++ b/python/ffsim/molecular_data.py
@@ -31,7 +31,6 @@
 import pyscf.mp
 import pyscf.symm
 import pyscf.tools
-from typing_extensions import deprecated
 
 from ffsim.hamiltonians import MolecularHamiltonian
 
@@ -194,31 +193,6 @@ def from_scf(
             hf_energy=hf_energy,
         )
 
-    @staticmethod
-    @deprecated(
-        "The from_mole method is deprecated. Instead, pass an SCF object directly to "
-        "from_scf."
-    )
-    def from_mole(
-        molecule: pyscf.gto.Mole,
-        active_space: Iterable[int] | None = None,
-        scf_func=pyscf.scf.RHF,
-    ) -> "MolecularData":
-        """Initialize a MolecularData object from a PySCF molecule.
-
-        .. warning::
-            This method is deprecated. Instead, pass an SCF object directly to
-            :func:`from_scf`.
-
-        Args:
-            molecule: The molecule.
-            active_space: An optional list of orbitals to use for the active space.
-            scf_func: The PySCF SCF function to use for the Hartree-Fock calculation.
-        """
-        hartree_fock = scf_func(molecule)
-        hartree_fock.run()
-        return MolecularData.from_scf(hartree_fock, active_space=active_space)
-
     def run_cisd(self, *, store_cisd_vec: bool = False) -> None:
         """Run CISD and store results."""
         cisd = pyscf.ci.CISD(self.scf.run())
diff --git a/python/ffsim/random/__init__.py b/python/ffsim/random/__init__.py
index 4c0457d1f..ba0c39dc5 100644
--- a/python/ffsim/random/__init__.py
+++ b/python/ffsim/random/__init__.py
@@ -22,7 +22,6 @@
     random_real_symmetric_matrix,
     random_special_orthogonal,
     random_state_vector,
-    random_statevector,
     random_t2_amplitudes,
     random_two_body_tensor,
     random_uccsd_restricted,
@@ -43,7 +42,6 @@
     "random_orthogonal",
     "random_real_symmetric_matrix",
     "random_special_orthogonal",
-    "random_statevector",
     "random_state_vector",
     "random_t2_amplitudes",
     "random_two_body_tensor",
diff --git a/python/ffsim/random/random.py b/python/ffsim/random/random.py
index d6ecf7055..8238bbb82 100644
--- a/python/ffsim/random/random.py
+++ b/python/ffsim/random/random.py
@@ -14,40 +14,11 @@
 from collections import defaultdict
 
 import numpy as np
-from typing_extensions import deprecated
 
 from ffsim import hamiltonians, operators, variational
 from ffsim.variational.util import validate_interaction_pairs
 
 
-@deprecated(
-    "ffsim.random.random_statevector is deprecated. "
-    "Instead, use ffsim.random.random_state_vector."
-)
-def random_statevector(dim: int, *, seed=None, dtype=complex) -> np.ndarray:
-    """Return a random state vector sampled from the uniform distribution.
-
-    .. warning::
-        This function is deprecated. Use :func:`ffsim.random.random_state_vector`
-        instead.
-
-    Args:
-        dim: The dimension of the state vector.
-        seed: A seed to initialize the pseudorandom number generator.
-            Should be a valid input to ``np.random.default_rng``.
-        dtype: The data type to use for the result.
-
-    Returns:
-        The sampled state vector.
-    """
-    rng = np.random.default_rng(seed)
-    vec = rng.standard_normal(dim).astype(dtype, copy=False)
-    if np.issubdtype(dtype, np.complexfloating):
-        vec += 1j * rng.standard_normal(dim).astype(dtype, copy=False)
-    vec /= np.linalg.norm(vec)
-    return vec
-
-
 def random_state_vector(dim: int, *, seed=None, dtype=complex) -> np.ndarray:
     """Return a random state vector sampled from the uniform distribution.
 
diff --git a/python/ffsim/states/__init__.py b/python/ffsim/states/__init__.py
index 8140d2211..c9bd1421a 100644
--- a/python/ffsim/states/__init__.py
+++ b/python/ffsim/states/__init__.py
@@ -13,25 +13,21 @@
 from ffsim.states.bitstring import (
     BitstringType,
     addresses_to_strings,
-    indices_to_strings,
     strings_to_addresses,
-    strings_to_indices,
 )
 from ffsim.states.product_state_sum import ProductStateSum
-from ffsim.states.rdm import rdm, rdms
+from ffsim.states.rdm import rdms
 from ffsim.states.sample_slater import sample_slater_determinant
 from ffsim.states.slater import (
     hartree_fock_state,
     slater_determinant,
     slater_determinant_amplitudes,
-    slater_determinant_rdm,
     slater_determinant_rdms,
 )
 from ffsim.states.states import (
     StateVector,
     dim,
     dims,
-    one_hot,
     sample_state_vector,
     spin_square,
 )
@@ -47,17 +43,12 @@
     "expectation_one_body_power",
     "expectation_one_body_product",
     "hartree_fock_state",
-    "indices_to_strings",
-    "one_hot",
-    "rdm",
     "rdms",
     "sample_slater_determinant",
     "sample_state_vector",
     "slater_determinant",
     "slater_determinant_amplitudes",
-    "slater_determinant_rdm",
     "slater_determinant_rdms",
     "spin_square",
     "strings_to_addresses",
-    "strings_to_indices",
 ]
diff --git a/python/ffsim/states/bitstring.py b/python/ffsim/states/bitstring.py
index cc2f2bc65..dc5b2acf3 100644
--- a/python/ffsim/states/bitstring.py
+++ b/python/ffsim/states/bitstring.py
@@ -19,7 +19,6 @@
 
 import numpy as np
 from pyscf.fci import cistring
-from typing_extensions import deprecated
 
 
 class BitstringType(Enum):
@@ -46,90 +45,6 @@ class BitstringType(Enum):
     """Bit array."""
 
 
-@deprecated(
-    "ffsim.indices_to_strings is deprecated. "
-    "Instead, use ffsim.addresses_to_strings."
-)
-def indices_to_strings(
-    indices: Sequence[int] | np.ndarray,
-    norb: int,
-    nelec: int | tuple[int, int],
-    concatenate: bool = True,
-    bitstring_type: BitstringType = BitstringType.STRING,
-):
-    """Convert state vector indices to bitstrings.
-
-    .. warning::
-        This function is deprecated. Use :class:`ffsim.addresses_to_strings` instead.
-
-    Example:
-
-    .. code::
-
-        import ffsim
-
-        norb = 3
-        nelec = (2, 1)
-        dim = ffsim.dim(norb, nelec)
-        ffsim.indices_to_strings(range(dim), norb, nelec)
-        # output:
-        # ['001011',
-        #  '010011',
-        #  '100011',
-        #  '001101',
-        #  '010101',
-        #  '100101',
-        #  '001110',
-        #  '010110',
-        #  '100110']
-
-    Args:
-        indices: The state vector indices to convert to bitstrings.
-        norb: The number of spatial orbitals.
-        nelec: Either a single integer representing the number of fermions for a
-            spinless system, or a pair of integers storing the numbers of spin alpha
-            and spin beta fermions.
-        bitstring_type: The desired type of bitstring output.
-        concatenate: Whether to concatenate the spin-alpha and spin-beta parts of the
-            bitstrings. If True, then a single list of concatenated bitstrings is
-            returned. The strings are concatenated in the order :math:`s_b s_a`,
-            that is, the alpha string appears on the right.
-            If False, then two lists are returned, ``(strings_a, strings_b)``. Note that
-            the list of alpha strings appears first, that is, on the left.
-            In the spinless case (when `nelec` is an integer), this argument is ignored.
-    """
-    if isinstance(nelec, int):
-        # Spinless case
-        return convert_bitstring_type(
-            list(cistring.addrs2str(norb=norb, nelec=nelec, addrs=indices)),
-            input_type=BitstringType.INT,
-            output_type=bitstring_type,
-            length=norb,
-        )
-
-    # Spinful case
-    n_alpha, n_beta = nelec
-    dim_b = math.comb(norb, n_beta)
-    indices_a, indices_b = np.divmod(indices, dim_b)
-    strings_a = convert_bitstring_type(
-        list(cistring.addrs2str(norb=norb, nelec=n_alpha, addrs=indices_a)),
-        input_type=BitstringType.INT,
-        output_type=bitstring_type,
-        length=norb,
-    )
-    strings_b = convert_bitstring_type(
-        list(cistring.addrs2str(norb=norb, nelec=n_beta, addrs=indices_b)),
-        input_type=BitstringType.INT,
-        output_type=bitstring_type,
-        length=norb,
-    )
-    if concatenate:
-        return concatenate_bitstrings(
-            strings_a, strings_b, bitstring_type=bitstring_type, length=norb
-        )
-    return strings_a, strings_b
-
-
 def convert_bitstring_type(
     strings: Sequence[str] | Sequence[int] | np.ndarray,
     input_type: BitstringType,
@@ -244,59 +159,6 @@ def concatenate_bitstrings(
         return np.concatenate([strings_b, strings_a], axis=1)
 
 
-@deprecated(
-    "ffsim.strings_to_indices is deprecated. "
-    "Instead, use ffsim.strings_to_addresses."
-)
-def strings_to_indices(
-    strings: Sequence[str], norb: int, nelec: int | tuple[int, int]
-) -> np.ndarray:
-    """Convert bitstrings to state vector indices.
-
-    .. warning::
-        This function is deprecated. Use :class:`ffsim.strings_to_addresses` instead.
-
-    Example:
-
-    .. code::
-
-        import ffsim
-
-        norb = 3
-        nelec = (2, 1)
-        dim = ffsim.dim(norb, nelec)
-        ffsim.strings_to_indices(
-            [
-                "001011",
-                "010011",
-                "100011",
-                "001101",
-                "010101",
-                "100101",
-                "001110",
-                "010110",
-                "100110",
-            ],
-            norb,
-            nelec,
-        )
-        # output:
-        # array([0, 1, 2, 3, 4, 5, 6, 7, 8], dtype=int32)
-    """
-    if isinstance(nelec, int):
-        return cistring.strs2addr(
-            norb=norb, nelec=nelec, strings=[int(s, base=2) for s in strings]
-        )
-
-    n_alpha, n_beta = nelec
-    strings_a = [int(s[norb:], base=2) for s in strings]
-    strings_b = [int(s[:norb], base=2) for s in strings]
-    addrs_a = cistring.strs2addr(norb=norb, nelec=n_alpha, strings=strings_a)
-    addrs_b = cistring.strs2addr(norb=norb, nelec=n_beta, strings=strings_b)
-    dim_b = math.comb(norb, n_beta)
-    return addrs_a * dim_b + addrs_b
-
-
 def addresses_to_strings(
     addresses: Sequence[int] | np.ndarray,
     norb: int,
diff --git a/python/ffsim/states/rdm.py b/python/ffsim/states/rdm.py
index 28754daab..b5c0043fe 100644
--- a/python/ffsim/states/rdm.py
+++ b/python/ffsim/states/rdm.py
@@ -13,7 +13,6 @@
 from __future__ import annotations
 
 import numpy as np
-import scipy.linalg
 from pyscf.fci.direct_spin1 import (
     make_rdm1,
     make_rdm1s,
@@ -24,210 +23,10 @@
     trans_rdm12,
     trans_rdm12s,
 )
-from typing_extensions import deprecated
 
 from ffsim.cistring import gen_linkstr_index
 
 
-@deprecated("ffsim.rdm is deprecated. Instead, use ffsim.rdms.")
-def rdm(
-    vec: np.ndarray,
-    norb: int,
-    nelec: tuple[int, int],
-    rank: int = 1,
-    spin_summed: bool = True,
-    reordered: bool = True,
-    return_lower_ranks: bool = True,
-) -> np.ndarray | tuple[np.ndarray, ...]:
-    """Return the reduced density matrix (RDM) or matrices of a state vector.
-
-    .. warning::
-        This function is deprecated. Use :func:`ffsim.rdms` instead.
-
-    The rank 1 RDM is defined as follows:
-
-    .. code::
-
-        rdm1[p, q] = ⟨p+ q⟩
-
-    The definition of higher-rank RDMs depends on the ``reordered`` argument, which
-    defaults to True.
-
-    **reordered = True**
-
-    The reordered RDMs are defined as follows:
-
-    .. code::
-
-        rdm2[p, q, r, s] = ⟨p+ r+ s q⟩
-        rdm3[p, q, r, s, t, u] = ⟨p+ r+ t+ u s q⟩
-        rdm4[p, q, r, s, t, u, v, w] = ⟨p+ r+ t+ v+ w u s q⟩
-
-    **reordered = False**
-
-    If reordered is set to False, the RDMs are defined as follows:
-
-    .. code::
-
-        rdm2[p, q, r, s] = ⟨p+ q r+ s⟩
-        rdm3[p, q, r, s, t, u] = ⟨p+ q r+ s t+ u⟩
-        rdm4[p, q, r, s, t, u, v, w] = ⟨p+ q r+ s t+ u v+ w⟩
-
-    Note:
-        Currently, only ranks 1 and 2 are supported.
-
-    Args:
-        vec: The state vector whose reduced density matrix is desired.
-        norb: The number of spatial orbitals.
-        nelec: The number of alpha and beta electrons.
-        rank: The rank of the reduced density matrix.
-        spin_summed: Whether to sum over the spin index.
-        reordered: Whether to reorder the indices of the reduced density matrix.
-        return_lower_ranks: Whether to return lower rank RDMs in addition to the
-            specified rank. If True, then this function returns all RDMs up to and
-            including the specified rank, in increasing order of rank. For example,
-            if `rank=2` then a tuple `(rdm1, rdm2)` is returned.
-
-    Returns:
-        The reduced density matrix or matrices. If `return_lower_ranks` is False,
-        then a single matrix is returned. If `return_lower_ranks` is True, then
-        a `rank`-length tuple of matrices is returned, containing the RDMs up to
-        the specified rank in increasing order of rank.
-    """
-    n_alpha, n_beta = nelec
-    link_index_a = gen_linkstr_index(range(norb), n_alpha)
-    link_index_b = gen_linkstr_index(range(norb), n_beta)
-    link_index = (link_index_a, link_index_b)
-    if rank == 1:
-        if spin_summed:
-            return _rdm1_spin_summed(vec, norb, nelec, link_index)
-        else:
-            return _rdm1(vec, norb, nelec, link_index)
-    if rank == 2:
-        if spin_summed:
-            return _rdm2_spin_summed_deprecated(
-                vec, norb, nelec, reordered, link_index, return_lower_ranks
-            )
-        else:
-            return _rdm2(vec, norb, nelec, reordered, link_index, return_lower_ranks)
-    raise NotImplementedError(
-        f"Computing the rank {rank} reduced density matrix is currently not supported."
-    )
-
-
-def _rdm2_spin_summed_deprecated(
-    vec: np.ndarray,
-    norb: int,
-    nelec: tuple[int, int],
-    reordered: bool,
-    link_index: tuple[np.ndarray, np.ndarray] | None,
-    return_lower_ranks: bool,
-) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
-    rdm1_real, rdm2_real = make_rdm12(
-        vec.real, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    rdm1_imag, rdm2_imag = make_rdm12(
-        vec.imag, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    trans_rdm1_real_imag, trans_rdm2_real_imag = trans_rdm12(
-        vec.real, vec.imag, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    rdm2 = _assemble_rdm2_spin_summed(rdm2_real, rdm2_imag, trans_rdm2_real_imag)
-    if not return_lower_ranks:
-        return rdm2
-    rdm1 = _assemble_rdm1_spin_summed(rdm1_real, rdm1_imag, trans_rdm1_real_imag)
-    return rdm1, rdm2
-
-
-def _rdm1(
-    vec: np.ndarray,
-    norb: int,
-    nelec: tuple[int, int],
-    link_index: tuple[np.ndarray, np.ndarray] | None,
-) -> np.ndarray:
-    rdms1_real = make_rdm1s(vec.real, norb, nelec, link_index=link_index)
-    rdms1_imag = make_rdm1s(vec.imag, norb, nelec, link_index=link_index)
-    trans_rdms1_real_imag = trans_rdm1s(
-        vec.real, vec.imag, norb, nelec, link_index=link_index
-    )
-    trans_rdms1_imag_real = trans_rdm1s(
-        vec.imag, vec.real, norb, nelec, link_index=link_index
-    )
-    return _assemble_rdm1(
-        rdms1_real, rdms1_imag, trans_rdms1_real_imag, trans_rdms1_imag_real
-    )
-
-
-def _assemble_rdm1(
-    rdms1_real: tuple[np.ndarray, np.ndarray],
-    rdms1_imag: tuple[np.ndarray, np.ndarray],
-    trans_rdms1_real_imag: tuple[np.ndarray, np.ndarray],
-    trans_rdms1_imag_real: tuple[np.ndarray, np.ndarray],
-) -> np.ndarray:
-    rdms1 = np.stack(rdms1_real).astype(complex)
-    rdms1 += np.stack(rdms1_imag)
-    # use minus sign for real-imag and plus sign for imag-real because
-    # the rdm1 convention in pyscf is transposed
-    rdms1 -= 1j * np.stack(trans_rdms1_real_imag)
-    rdms1 += 1j * np.stack(trans_rdms1_imag_real)
-    return scipy.linalg.block_diag(*rdms1)
-
-
-def _rdm2(
-    vec: np.ndarray,
-    norb: int,
-    nelec: tuple[int, int],
-    reordered: bool,
-    link_index: tuple[np.ndarray, np.ndarray] | None,
-    return_lower_ranks: bool,
-) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
-    rdms1_real, rdms2_real = make_rdm12s(
-        vec.real, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    rdms1_imag, rdms2_imag = make_rdm12s(
-        vec.imag, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    trans_rdms1_real_imag, trans_rdms2_real_imag = trans_rdm12s(
-        vec.real, vec.imag, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    trans_rdms1_imag_real, trans_rdms2_imag_real = trans_rdm12s(
-        vec.imag, vec.real, norb, nelec, reorder=reordered, link_index=link_index
-    )
-    rdm2 = _assemble_rdm2(
-        rdms2_real, rdms2_imag, trans_rdms2_real_imag, trans_rdms2_imag_real
-    )
-    if not return_lower_ranks:
-        return rdm2
-    rdm1 = _assemble_rdm1(
-        rdms1_real, rdms1_imag, trans_rdms1_real_imag, trans_rdms1_imag_real
-    )
-    return rdm1, rdm2
-
-
-def _assemble_rdm2(
-    rdms2_real: tuple[np.ndarray, np.ndarray, np.ndarray],
-    rdms2_imag: tuple[np.ndarray, np.ndarray, np.ndarray],
-    trans_rdms2_real_imag: tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
-    trans_rdms2_imag_real: tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
-) -> np.ndarray:
-    rdms2 = np.stack(rdms2_real).astype(complex)
-    rdms2 += np.stack(rdms2_imag)
-    # rdms2 is currently [rdm_aa, rdm_ab, rdm_bb]
-    # the following line transforms it into [rdm_aa, rdm_ab, rdm_ba, rdm_bb]
-    rdms2 = np.insert(rdms2, 2, rdms2[1].transpose(2, 3, 0, 1), axis=0)
-    rdms2 += 1j * np.stack(trans_rdms2_real_imag)
-    rdms2 -= 1j * np.stack(trans_rdms2_imag_real)
-
-    norb, _, _, _ = rdms2_real[0].shape
-    rdm2 = np.zeros((2 * norb, 2 * norb, 2 * norb, 2 * norb), dtype=complex)
-    rdm_aa, rdm_ab, rdm_ba, rdm_bb = rdms2
-    rdm2[:norb, :norb, :norb, :norb] = rdm_aa
-    rdm2[:norb, :norb, norb:, norb:] = rdm_ab
-    rdm2[norb:, norb:, :norb, :norb] = rdm_ba
-    rdm2[norb:, norb:, norb:, norb:] = rdm_bb
-    return rdm2
-
-
 def rdms(
     vec: np.ndarray,
     norb: int,
@@ -237,7 +36,7 @@ def rdms(
     spin_summed: bool = False,
     reorder: bool = True,
 ) -> np.ndarray | tuple[np.ndarray, ...]:
-    """Return the spin-separated reduced density matrices of a state vector.
+    """Return the reduced density matrices of a state vector.
 
     The rank 1 RDM is defined as follows:
 
diff --git a/python/ffsim/states/slater.py b/python/ffsim/states/slater.py
index 77ec95672..feb290d08 100644
--- a/python/ffsim/states/slater.py
+++ b/python/ffsim/states/slater.py
@@ -18,7 +18,6 @@
 import numpy as np
 import scipy.linalg
 from pyscf.fci import cistring
-from typing_extensions import deprecated
 
 from ffsim import linalg
 from ffsim.gates.orbital_rotation import apply_orbital_rotation
@@ -128,76 +127,6 @@ def hartree_fock_state(norb: int, nelec: int | tuple[int, int]) -> np.ndarray:
     return slater_determinant(norb, occupied_orbitals=(range(n_alpha), range(n_beta)))
 
 
-@deprecated(
-    "ffsim.slater_determinant_rdm is deprecated. "
-    "Instead, use ffsim.slater_determinant_rdms."
-)
-def slater_determinant_rdm(
-    norb: int,
-    occupied_orbitals: tuple[Sequence[int], Sequence[int]],
-    orbital_rotation: np.ndarray
-    | tuple[np.ndarray | None, np.ndarray | None]
-    | None = None,
-    rank: int = 1,
-    spin_summed: bool = True,
-) -> np.ndarray:
-    """Return the reduced density matrix of a `Slater determinant`_.
-
-    .. warning::
-        This function is deprecated. Use :func:`ffsim.slater_determinant_rdms` instead.
-
-    Note:
-        Currently, only rank 1 is supported.
-
-    Args:
-        norb: The number of spatial orbitals.
-        occupied_orbitals: The occupied orbitals in the electronic configuration.
-            This is a pair of lists of integers, where the first list specifies the
-            spin alpha orbitals and the second list specifies the spin beta
-            orbitals.
-        orbital_rotation: The optional orbital rotation.
-            You can pass either a single Numpy array specifying the orbital rotation
-            to apply to both spin sectors, or you can pass a pair of Numpy arrays
-            specifying independent orbital rotations for spin alpha and spin beta.
-            If passing a pair, you can use ``None`` for one of the
-            values in the pair to indicate that no operation should be applied to that
-            spin sector.
-        rank: The rank of the reduced density matrix. I.e., rank 1 corresponds to the
-            one-particle RDM, rank 2 corresponds to the 2-particle RDM, etc.
-        spin_summed: Whether to sum over the spin index.
-
-    Returns:
-        The reduced density matrix of the Slater determinant.
-
-    .. _Slater determinant: ffsim.html#ffsim.slater_determinant
-    """
-    if rank == 1:
-        rdm_a = np.zeros((norb, norb), dtype=complex)
-        rdm_b = np.zeros((norb, norb), dtype=complex)
-        alpha_orbitals = np.array(occupied_orbitals[0])
-        beta_orbitals = np.array(occupied_orbitals[1])
-        if len(alpha_orbitals):
-            rdm_a[(alpha_orbitals, alpha_orbitals)] = 1
-        if len(beta_orbitals):
-            rdm_b[(beta_orbitals, beta_orbitals)] = 1
-        if orbital_rotation is not None:
-            if isinstance(orbital_rotation, np.ndarray):
-                orbital_rotation_a: np.ndarray | None = orbital_rotation
-                orbital_rotation_b: np.ndarray | None = orbital_rotation
-            else:
-                orbital_rotation_a, orbital_rotation_b = orbital_rotation
-            if orbital_rotation_a is not None:
-                rdm_a = orbital_rotation_a.conj() @ rdm_a @ orbital_rotation_a.T
-            if orbital_rotation_b is not None:
-                rdm_b = orbital_rotation_b.conj() @ rdm_b @ orbital_rotation_b.T
-        if spin_summed:
-            return rdm_a + rdm_b
-        return scipy.linalg.block_diag(rdm_a, rdm_b)
-    raise NotImplementedError(
-        f"Returning the rank {rank} reduced density matrix is currently not supported."
-    )
-
-
 @overload
 def slater_determinant_rdms(
     norb: int,
diff --git a/python/ffsim/states/states.py b/python/ffsim/states/states.py
index dd2faf89d..4bdadfd3f 100644
--- a/python/ffsim/states/states.py
+++ b/python/ffsim/states/states.py
@@ -19,7 +19,6 @@
 
 import numpy as np
 from pyscf.fci.spin_op import contract_ss
-from typing_extensions import deprecated
 
 from ffsim.states.bitstring import (
     BitstringType,
@@ -93,28 +92,6 @@ def dim(norb: int, nelec: int | tuple[int, int]) -> int:
     return math.comb(norb, n_alpha) * math.comb(norb, n_beta)
 
 
-@deprecated(
-    "Using one_hot from the ffsim namespace is deprecated. "
-    "Instead, use ffsim.linalg.one_hot."
-)
-def one_hot(shape: int | tuple[int, ...], index, *, dtype=complex):
-    """Return an array of all zeros except for a one at a specified index.
-
-    .. warning::
-        This function is deprecated. Use :func:`ffsim.linalg.one_hot` instead.
-
-    Args:
-        shape: The desired shape of the array.
-        index: The index at which to place a one.
-
-    Returns:
-        The one-hot vector.
-    """
-    vec = np.zeros(shape, dtype=dtype)
-    vec[index] = 1
-    return vec
-
-
 # source: pyscf.fci.spin_op.spin_square0
 # modified to support complex wavefunction
 def spin_square(fcivec: np.ndarray, norb: int, nelec: tuple[int, int]):
diff --git a/tests/python/hamiltonians/molecular_hamiltonian_test.py b/tests/python/hamiltonians/molecular_hamiltonian_test.py
index 3b594970a..58524d557 100644
--- a/tests/python/hamiltonians/molecular_hamiltonian_test.py
+++ b/tests/python/hamiltonians/molecular_hamiltonian_test.py
@@ -12,8 +12,6 @@
 
 from __future__ import annotations
 
-import pathlib
-
 import numpy as np
 import pyscf
 import pyscf.mcscf
@@ -164,28 +162,3 @@ def test_rotated():
     original_expectation = np.vdot(vec, linop @ vec)
     rotated_expectation = np.vdot(rotated_vec, linop_rotated @ rotated_vec)
     np.testing.assert_allclose(original_expectation, rotated_expectation)
-
-
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_from_fcidump(tmp_path: pathlib.Path):
-    """Test loading from FCIDUMP."""
-    mol = pyscf.gto.Mole()
-    mol.build(
-        atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
-        basis="sto-6g",
-        symmetry="Dooh",
-    )
-    n_frozen = pyscf.data.elements.chemcore(mol)
-    active_space = range(n_frozen, mol.nao_nr())
-    scf = pyscf.scf.RHF(mol).run()
-    mol_data = ffsim.MolecularData.from_scf(scf, active_space=active_space)
-    pyscf.tools.fcidump.from_integrals(
-        tmp_path / "test.fcidump",
-        h1e=mol_data.one_body_integrals,
-        h2e=mol_data.two_body_integrals,
-        nuc=mol_data.core_energy,
-        nmo=mol_data.norb,
-        nelec=mol_data.nelec,
-    )
-    mol_ham = ffsim.MolecularHamiltonian.from_fcidump(tmp_path / "test.fcidump")
-    assert ffsim.approx_eq(mol_ham, mol_data.hamiltonian)
diff --git a/tests/python/qiskit/sampler_test.py b/tests/python/qiskit/sampler_test.py
index 33fe17515..1257c3399 100644
--- a/tests/python/qiskit/sampler_test.py
+++ b/tests/python/qiskit/sampler_test.py
@@ -302,8 +302,6 @@ def test_global_depolarizing(
     assert np.allclose(fidelity, expected_fidelity, rtol=1e-2, atol=1e-3)
 
 
-# TODO remove after removing UCJOperatorJW
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
 def test_reproducible_with_seed():
     """Test sampler with random gates."""
     rng = np.random.default_rng(14062)
diff --git a/tests/python/states/bitstring_test.py b/tests/python/states/bitstring_test.py
index 47d75e695..4cad6562e 100644
--- a/tests/python/states/bitstring_test.py
+++ b/tests/python/states/bitstring_test.py
@@ -19,242 +19,6 @@
 from ffsim.states.bitstring import convert_bitstring_type
 
 
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_indices_to_strings_string():
-    """Test converting statevector indices to strings, output type string."""
-    norb = 3
-    nelec: int | tuple[int, int] = 2
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(range(dim), norb, nelec)
-    assert strings == [
-        "011",
-        "101",
-        "110",
-    ]
-    strings = ffsim.indices_to_strings(range(dim), norb, nelec, concatenate=False)
-    assert strings == [
-        "011",
-        "101",
-        "110",
-    ]
-
-    norb = 3
-    nelec = (2, 1)
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(range(dim), norb, nelec)
-    assert strings == [
-        "001011",
-        "010011",
-        "100011",
-        "001101",
-        "010101",
-        "100101",
-        "001110",
-        "010110",
-        "100110",
-    ]
-    strings_a, strings_b = ffsim.indices_to_strings(
-        range(dim), norb, nelec, concatenate=False
-    )
-    assert strings_a == [
-        "011",
-        "011",
-        "011",
-        "101",
-        "101",
-        "101",
-        "110",
-        "110",
-        "110",
-    ]
-    assert strings_b == [
-        "001",
-        "010",
-        "100",
-        "001",
-        "010",
-        "100",
-        "001",
-        "010",
-        "100",
-    ]
-
-
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_indices_to_strings_int():
-    """Test converting statevector indices to strings, output type int."""
-    norb = 3
-    nelec: int | tuple[int, int] = 2
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(
-        range(dim), norb, nelec, bitstring_type=ffsim.BitstringType.INT
-    )
-    assert strings == [
-        0b011,
-        0b101,
-        0b110,
-    ]
-    strings = ffsim.indices_to_strings(
-        range(dim),
-        norb,
-        nelec,
-        concatenate=False,
-        bitstring_type=ffsim.BitstringType.INT,
-    )
-    assert strings == [
-        0b011,
-        0b101,
-        0b110,
-    ]
-
-    norb = 3
-    nelec = (2, 1)
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(
-        range(dim), norb, nelec, bitstring_type=ffsim.BitstringType.INT
-    )
-    assert strings == [
-        0b001011,
-        0b010011,
-        0b100011,
-        0b001101,
-        0b010101,
-        0b100101,
-        0b001110,
-        0b010110,
-        0b100110,
-    ]
-    strings_a, strings_b = ffsim.indices_to_strings(
-        range(dim),
-        norb,
-        nelec,
-        concatenate=False,
-        bitstring_type=ffsim.BitstringType.INT,
-    )
-    assert strings_a == [
-        0b011,
-        0b011,
-        0b011,
-        0b101,
-        0b101,
-        0b101,
-        0b110,
-        0b110,
-        0b110,
-    ]
-    assert strings_b == [
-        0b001,
-        0b010,
-        0b100,
-        0b001,
-        0b010,
-        0b100,
-        0b001,
-        0b010,
-        0b100,
-    ]
-
-
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_indices_to_strings_bit_array():
-    """Test converting statevector indices to strings, output type bit array."""
-    norb = 3
-    nelec: int | tuple[int, int] = 2
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(
-        range(dim), norb, nelec, bitstring_type=ffsim.BitstringType.BIT_ARRAY
-    )
-    np.testing.assert_array_equal(
-        strings,
-        np.array(
-            [
-                [False, True, True],
-                [True, False, True],
-                [True, True, False],
-            ]
-        ),
-    )
-    strings = ffsim.indices_to_strings(
-        range(dim),
-        norb,
-        nelec,
-        concatenate=False,
-        bitstring_type=ffsim.BitstringType.BIT_ARRAY,
-    )
-    np.testing.assert_array_equal(
-        strings,
-        np.array(
-            [
-                [False, True, True],
-                [True, False, True],
-                [True, True, False],
-            ]
-        ),
-    )
-
-    norb = 3
-    nelec = (2, 1)
-    dim = ffsim.dim(norb, nelec)
-    strings = ffsim.indices_to_strings(
-        range(dim), norb, nelec, bitstring_type=ffsim.BitstringType.BIT_ARRAY
-    )
-    np.testing.assert_array_equal(
-        strings,
-        np.array(
-            [
-                [False, False, True, False, True, True],
-                [False, True, False, False, True, True],
-                [True, False, False, False, True, True],
-                [False, False, True, True, False, True],
-                [False, True, False, True, False, True],
-                [True, False, False, True, False, True],
-                [False, False, True, True, True, False],
-                [False, True, False, True, True, False],
-                [True, False, False, True, True, False],
-            ]
-        ),
-    )
-    strings_a, strings_b = ffsim.indices_to_strings(
-        range(dim),
-        norb,
-        nelec,
-        concatenate=False,
-        bitstring_type=ffsim.BitstringType.BIT_ARRAY,
-    )
-    np.testing.assert_array_equal(
-        strings_a,
-        np.array(
-            [
-                [False, True, True],
-                [False, True, True],
-                [False, True, True],
-                [True, False, True],
-                [True, False, True],
-                [True, False, True],
-                [True, True, False],
-                [True, True, False],
-                [True, True, False],
-            ]
-        ),
-    )
-    np.testing.assert_array_equal(
-        strings_b,
-        np.array(
-            [
-                [False, False, True],
-                [False, True, False],
-                [True, False, False],
-                [False, False, True],
-                [False, True, False],
-                [True, False, False],
-                [False, False, True],
-                [False, True, False],
-                [True, False, False],
-            ]
-        ),
-    )
-
-
 def test_addresses_to_strings_int_spinless():
     """Test converting statevector addresses to strings, int output, spinless."""
     norb = 3
diff --git a/tests/python/states/rdm_test.py b/tests/python/states/rdm_test.py
index 7d2314b8a..dd30b4e23 100644
--- a/tests/python/states/rdm_test.py
+++ b/tests/python/states/rdm_test.py
@@ -16,7 +16,6 @@
 
 import numpy as np
 import pytest
-import scipy.linalg
 
 import ffsim
 
@@ -208,142 +207,3 @@ def _rdm2_spin_summed_reordered(
     """Compute spin-summed reordered 2-RDM directly from its definition."""
     rdm_aa, rdm_ab, rdm_bb = _rdm2s_reordered(vec, norb, nelec)
     return rdm_aa + rdm_ab + rdm_ab.transpose(2, 3, 0, 1) + rdm_bb
-
-
-@pytest.mark.parametrize(
-    "norb, nelec, spin_summed",
-    [
-        (4, (2, 2), True),
-        (4, (2, 2), False),
-        (3, (1, 2), True),
-        (3, (1, 2), False),
-        (2, (0, 1), True),
-        (2, (0, 1), False),
-        (1, (0, 0), True),
-        (1, (0, 0), False),
-    ],
-)
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_rdm_1(norb: int, nelec: tuple[int, int], spin_summed: bool):
-    """Test computing 1-RDM."""
-    func = {
-        # spin_summed: function
-        False: _rdm1,
-        True: _rdm1_spin_summed,
-    }
-    rng = np.random.default_rng()
-    vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
-
-    rdm = ffsim.rdm(vec, norb, nelec, spin_summed=spin_summed)
-    expected = func[spin_summed](vec, norb, nelec)
-    np.testing.assert_allclose(rdm, expected, atol=1e-12)
-
-
-@pytest.mark.parametrize(
-    "norb, nelec, spin_summed, reordered",
-    [
-        (4, (2, 2), True, True),
-        (4, (2, 2), True, False),
-        (4, (2, 2), False, True),
-        (4, (2, 2), False, False),
-        (3, (1, 2), True, True),
-        (3, (1, 2), True, False),
-        (3, (1, 2), False, True),
-        (3, (1, 2), False, False),
-        (2, (0, 1), True, True),
-        (2, (0, 1), True, False),
-        (2, (0, 1), False, True),
-        (2, (0, 1), False, False),
-        (1, (0, 0), True, True),
-        (1, (0, 0), True, False),
-        (1, (0, 0), False, True),
-        (1, (0, 0), False, False),
-    ],
-)
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_rdm_2(norb: int, nelec: tuple[int, int], spin_summed: bool, reordered: bool):
-    """Test computing 1- and 2-RDMs."""
-    func = {
-        # (spin_summed, reorder): (rdm1_function, rdm2_function)
-        (False, False): (_rdm1, _rdm2),
-        (False, True): (_rdm1, _rdm2_reordered),
-        (True, False): (_rdm1_spin_summed, _rdm2_spin_summed),
-        (True, True): (_rdm1_spin_summed, _rdm2_spin_summed_reordered),
-    }
-    rng = np.random.default_rng()
-    vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
-
-    rdm1, rdm2 = ffsim.rdm(
-        vec,
-        norb,
-        nelec,
-        rank=2,
-        spin_summed=spin_summed,
-        reordered=reordered,
-    )
-    rdm1_func, rdm2_func = func[spin_summed, reordered]
-    expected_rdm1 = rdm1_func(vec, norb, nelec)
-    expected_rdm2 = rdm2_func(vec, norb, nelec)
-    np.testing.assert_allclose(rdm1, expected_rdm1, atol=1e-12)
-    np.testing.assert_allclose(rdm2, expected_rdm2, atol=1e-12)
-
-
-@pytest.mark.parametrize(
-    "norb, nelec, spin_summed, reordered, rank",
-    [
-        (4, (2, 2), True, True, 2),
-    ],
-)
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_no_lower_ranks(
-    norb: int, nelec: tuple[int, int], spin_summed: bool, reordered: bool, rank: int
-):
-    """Test computing higher-rank RDM without returning lower ranks."""
-    func = {
-        # (spin_summed, reorder): function
-        (False, False, 2): _rdm2,
-        (False, True, 2): _rdm2_reordered,
-        (True, False, 2): _rdm2_spin_summed,
-        (True, True, 2): _rdm2_spin_summed_reordered,
-    }
-    rng = np.random.default_rng()
-    vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
-
-    rdm = ffsim.rdm(
-        vec,
-        norb,
-        nelec,
-        rank=rank,
-        spin_summed=spin_summed,
-        reordered=reordered,
-        return_lower_ranks=False,
-    )
-    expected = func[spin_summed, reordered, rank](vec, norb, nelec)
-    np.testing.assert_allclose(rdm, expected, atol=1e-12)
-
-
-def _rdm1(vec: np.ndarray, norb: int, nelec: tuple[int, int]) -> np.ndarray:
-    """Compute 1-RDM directly from its definition."""
-    return scipy.linalg.block_diag(*_rdm1s(vec, norb, nelec))
-
-
-def _rdm2(vec: np.ndarray, norb: int, nelec: tuple[int, int]) -> np.ndarray:
-    """Compute 2-RDM directly from its definition."""
-    rdm = np.zeros((2 * norb, 2 * norb, 2 * norb, 2 * norb), dtype=complex)
-    aa, ab, bb = _rdm2s(vec, norb, nelec)
-    rdm[:norb, :norb, :norb, :norb] = aa
-    rdm[:norb, :norb, norb:, norb:] = ab
-    rdm[norb:, norb:, :norb, :norb] = ab.transpose(2, 3, 0, 1)
-    rdm[norb:, norb:, norb:, norb:] = bb
-    return rdm
-
-
-def _rdm2_reordered(vec: np.ndarray, norb: int, nelec: tuple[int, int]) -> np.ndarray:
-    """Compute reordered 2-RDM directly from its definition."""
-    rdm = np.zeros((2 * norb, 2 * norb, 2 * norb, 2 * norb), dtype=complex)
-    aa, ab, bb = _rdm2s_reordered(vec, norb, nelec)
-    rdm[:norb, :norb, :norb, :norb] = aa
-    rdm[:norb, :norb, norb:, norb:] = ab
-    rdm[norb:, norb:, :norb, :norb] = ab.transpose(2, 3, 0, 1)
-    rdm[norb:, norb:, norb:, norb:] = bb
-    return rdm
diff --git a/tests/python/states/states_test.py b/tests/python/states/states_test.py
index c75f9b902..9e89af3ef 100644
--- a/tests/python/states/states_test.py
+++ b/tests/python/states/states_test.py
@@ -12,8 +12,6 @@
 
 from __future__ import annotations
 
-import itertools
-
 import numpy as np
 import pytest
 
@@ -182,39 +180,6 @@ def test_slater_determinant_rdm1s_spinless(norb: int, nelec: int):
     np.testing.assert_allclose(rdm, expected, atol=1e-12)
 
 
-@pytest.mark.parametrize(
-    "norb, nelec, spin_summed",
-    [
-        (norb, nelec, spin_summed)
-        for (norb, nelec), spin_summed in itertools.product(
-            ffsim.testing.generate_norb_nelec(range(5)), [False, True]
-        )
-    ],
-)
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_slater_determinant_one_rdm_same_rotation(
-    norb: int, nelec: tuple[int, int], spin_summed: bool
-):
-    """Test Slater determinant 1-RDM."""
-    rng = np.random.default_rng()
-
-    occupied_orbitals = ffsim.testing.random_occupied_orbitals(norb, nelec, seed=rng)
-    orbital_rotation = ffsim.random.random_unitary(norb, seed=rng)
-
-    vec = ffsim.slater_determinant(
-        norb, occupied_orbitals, orbital_rotation=orbital_rotation
-    )
-    rdm = ffsim.slater_determinant_rdm(
-        norb,
-        occupied_orbitals,
-        orbital_rotation=orbital_rotation,
-        spin_summed=spin_summed,
-    )
-    expected = ffsim.rdm(vec, norb, nelec, spin_summed=spin_summed)
-
-    np.testing.assert_allclose(rdm, expected, atol=1e-12)
-
-
 def test_sample_state_vector_spinful_string():
     """Test sampling state vector, spinful, output type string."""
     norb = 5
@@ -575,42 +540,6 @@ def test_sample_state_vector_spinless_bit_array():
     )
 
 
-@pytest.mark.parametrize(
-    "norb, nelec, spin_summed",
-    [
-        (norb, nelec, spin_summed)
-        for (norb, nelec), spin_summed in itertools.product(
-            ffsim.testing.generate_norb_nelec(range(5)), [False, True]
-        )
-    ],
-)
-@pytest.mark.filterwarnings("ignore::DeprecationWarning")
-def test_slater_determinant_one_rdm_diff_rotation(
-    norb: int, nelec: tuple[int, int], spin_summed: bool
-):
-    """Test Slater determinant 1-RDM."""
-    rng = np.random.default_rng()
-
-    occupied_orbitals = ffsim.testing.random_occupied_orbitals(norb, nelec, seed=rng)
-    orbital_rotation_a = ffsim.random.random_unitary(norb, seed=rng)
-    orbital_rotation_b = ffsim.random.random_unitary(norb, seed=rng)
-
-    vec = ffsim.slater_determinant(
-        norb,
-        occupied_orbitals,
-        orbital_rotation=(orbital_rotation_a, orbital_rotation_b),
-    )
-    rdm = ffsim.slater_determinant_rdm(
-        norb,
-        occupied_orbitals,
-        orbital_rotation=(orbital_rotation_a, orbital_rotation_b),
-        spin_summed=spin_summed,
-    )
-    expected = ffsim.rdm(vec, norb, nelec, spin_summed=spin_summed)
-
-    np.testing.assert_allclose(rdm, expected, atol=1e-12)
-
-
 def test_state_vector_array():
     """Test StateVector's __array__ method."""
     norb = 5