def kernel(self, eris: Optional[ERIsInputType] = None) -> float:
    """Run the coupled cluster calculation.

    Args:
        eris: Electron repulsion integrals.

    Returns:
        Correlation energy.
    """
    timer = util.Timer()

    # Get the ERIs:
    eris = self.get_eris(eris)

    # Get the amplitude guesses:
    amplitudes = self.amplitudes
    if not amplitudes:
        amplitudes = self.init_amps(eris=eris)

    # Get the initial energy:
    e_cc = self.energy(amplitudes=amplitudes, eris=eris)

    self.log.output("Solving for excitation amplitudes.")
    self.log.debug("")
    self.log.info(
        f"{ANSI.B}{'Iter':>4s} {'Energy (corr.)':>16s} {'Energy (tot.)':>18s} "
        f"{'Δ(Energy)':>13s} {'Δ(Ampl.)':>13s}{ANSI.R}"
    )
    self.log.info(f"{0:4d} {e_cc:16.10f} {e_cc + self.mf.e_tot:18.10f}")

    if not self.ansatz.is_one_shot:
        # Set up DIIS:
        diis = DIIS()
        diis.space = self.options.diis_space
        diis.damping = self.options.damping

        converged = False
        for niter in range(1, self.options.max_iter + 1):
            # Update the amplitudes, extrapolate with DIIS and calculate change:
            amplitudes_prev = amplitudes
            amplitudes = self.update_amps(amplitudes=amplitudes, eris=eris)
            vector = self.amplitudes_to_vector(amplitudes)
            vector = diis.update(vector)
            amplitudes = self.vector_to_amplitudes(vector)
            dt = np.max(np.abs(vector - self.amplitudes_to_vector(amplitudes_prev)))

            # Update the energy and calculate change:
            e_prev = e_cc
            e_cc = self.energy(amplitudes=amplitudes, eris=eris)
            de = abs(e_prev - e_cc)

            # Log the iteration:
            converged_e = bool(de < self.options.e_tol)
            converged_t = bool(dt < self.options.t_tol)
            self.log.info(
                f"{niter:4d} {e_cc:16.10f} {e_cc + self.mf.e_tot:18.10f}"
                f" {[ANSI.r, ANSI.g][bool(converged_e)]}{de:13.3e}{ANSI.R}"
                f" {[ANSI.r, ANSI.g][bool(converged_t)]}{dt:13.3e}{ANSI.R}"
            )

            # Check for convergence:
            converged = converged_e and converged_t
            if converged:
                self.log.debug("")
                self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
                break
        else:
            self.log.debug("")
            self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}")

        # Include perturbative correction if required:
        if self.ansatz.has_perturbative_correction:
            self.log.debug("")
            self.log.info("Computing perturbative energy correction.")
            e_pert = self.energy_perturbative(amplitudes=amplitudes, eris=eris)
            e_cc += e_pert
            self.log.info(f"E(pert) = {e_pert:.10f}")

    else:
        converged = True

    # Update attributes:
    self.e_corr = e_cc
    self.amplitudes = amplitudes
    self.converged = converged

    self.log.debug("")
    self.log.output(f"E(corr) = {self.e_corr:.10f}")
    self.log.output(f"E(tot)  = {self.e_tot:.10f}")
    self.log.debug("")
    self.log.debug("Time elapsed: %s", timer.format_time(timer()))
    self.log.debug("")

    return e_cc

def solve_lambda(
    self,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    eris: Optional[ERIsInputType] = None,
) -> None:
    """Solve for the lambda amplitudes.

    Args:
        amplitudes: Cluster amplitudes.
        eris: Electron repulsion integrals.
    """
    timer = util.Timer()

    # Get the ERIs:
    eris = self.get_eris(eris)

    # Get the amplitudes:
    amplitudes = self.amplitudes
    if not amplitudes:
        amplitudes = self.init_amps(eris=eris)

    # If needed, get the perturbative part of the lambda amplitudes:
    lambdas_pert = None
    if self.ansatz.has_perturbative_correction:
        lambdas_pert = self.update_lams(eris=eris, amplitudes=amplitudes, perturbative=True)

    # Get the initial lambda amplitudes:
    lambdas = self.lambdas
    if not lambdas:
        lambdas = self.init_lams(amplitudes=amplitudes)

    self.log.output("Solving for de-excitation (lambda) amplitudes.")
    self.log.debug("")
    self.log.info(f"{ANSI.B}{'Iter':>4s} {'Δ(Ampl.)':>13s}{ANSI.R}")

    # Set up DIIS:
    diis = DIIS()
    diis.space = self.options.diis_space
    diis.damping = self.options.damping

    converged = False
    for niter in range(1, self.options.max_iter + 1):
        # Update the lambda amplitudes, extrapolate with DIIS and calculate change:
        lambdas_prev = lambdas
        lambdas = self.update_lams(
            amplitudes=amplitudes,
            lambdas=lambdas,
            lambdas_pert=lambdas_pert,
            eris=eris,
        )
        vector = self.lambdas_to_vector(lambdas)
        vector = diis.update(vector)
        lambdas = self.vector_to_lambdas(vector)
        dl = np.max(np.abs(vector - self.lambdas_to_vector(lambdas_prev)))

        # Log the iteration:
        converged = bool(dl < self.options.t_tol)
        self.log.info(f"{niter:4d} {[ANSI.r, ANSI.g][converged]}{dl:13.3e}{ANSI.R}")

        # Check for convergence:
        if converged:
            self.log.debug("")
            self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
            break
    else:
        self.log.debug("")
        self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}")

    self.log.debug("")
    self.log.debug("Time elapsed: %s", timer.format_time(timer()))
    self.log.debug("")
    self.log.debug("")

    # Update attributes:
    self.lambdas = lambdas
    self.converged_lambda = converged

@abstractmethod
def ip_eom(self, **kwargs: Any) -> Any:
    """Get the IP-EOM object.

    Args:
        **kwargs: Additional keyword arguments.

    Returns:
        IP-EOM object.
    """
    pass

@abstractmethod
def ea_eom(self, **kwargs: Any) -> Any:
    """Get the EA-EOM object.

    Args:
        **kwargs: Additional keyword arguments.

    Returns:
        EA-EOM object.
    """
    pass

@abstractmethod
def ee_eom(self, **kwargs: Any) -> Any:
    """Get the EE-EOM object.

    Args:
        **kwargs: Additional keyword arguments.

    Returns:
        EE-EOM object.
    """
    pass

def brueckner(self, *args: Any, **kwargs: Any) -> float:
    """Run a Brueckner orbital coupled cluster calculation.

    The coupled cluster object will be update in-place.

    Args:
        *args: Arguments to pass to the Brueckner object.
        **kwargs: Keyword arguments to pass to the Brueckner object.

    Returns:
        Correlation energy.
    """
    bcc = self.Brueckner(self, *args, **kwargs)
    return bcc.kernel()

def write(self, file: str) -> None:
    """Write the EBCC object to a file.

    Args:
        file: File to write the object to.
    """
    writer = Dump(file)
    writer.write(self)

@classmethod
def read(cls, file: str, log: Optional[Logger] = None) -> BaseEBCC:
    """Read the EBCC object from a file.

    Args:
        file: File to read the object from.
        log: Logger to use for new object.

    Returns:
        EBCC object.
    """
    reader = Dump(file)
    return reader.read(cls=cls, log=log)

@abstractmethod
def init_amps(self, eris: Optional[ERIsInputType] = None) -> Namespace[SpinArrayType]:
    """Initialise the cluster amplitudes.

    Args:
        eris: Electron repulsion integrals.

    Returns:
        Initial cluster amplitudes.
    """
    pass

@abstractmethod
def init_lams(
    self, amplitudes: Optional[Namespace[SpinArrayType]] = None
) -> Namespace[SpinArrayType]:
    """Initialise the cluster lambda amplitudes.

    Args:
        amplitudes: Cluster amplitudes.

    Returns:
        Initial cluster lambda amplitudes.
    """
    pass

def energy(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
) -> float:
    """Calculate the correlation energy.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.

    Returns:
        Correlation energy.
    """
    func, kwargs = self._load_function(
        "energy",
        eris=eris,
        amplitudes=amplitudes,
    )
    res: float = func(**kwargs).real.item()
    return astype(res, float)

def energy_perturbative(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
) -> float:
    """Calculate the perturbative correction to the correlation energy.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.

    Returns:
        Perturbative energy correction.
    """
    func, kwargs = self._load_function(
        "energy_perturbative",
        eris=eris,
        amplitudes=amplitudes,
        lambdas=lambdas,
    )
    res: float = func(**kwargs).real.item()
    return astype(res, float)

@abstractmethod
def update_amps(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
) -> Namespace[SpinArrayType]:
    """Update the cluster amplitudes.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.

    Returns:
        Updated cluster amplitudes.
    """
    pass

@abstractmethod
def update_lams(
    self,
    eris: ERIsInputType = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
    lambdas_pert: Optional[Namespace[SpinArrayType]] = None,
    perturbative: bool = False,
) -> Namespace[SpinArrayType]:
    """Update the cluster lambda amplitudes.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.
        lambdas_pert: Perturbative cluster lambda amplitudes.
        perturbative: Flag to include perturbative correction.

    Returns:
        Updated cluster lambda amplitudes.
    """
    pass

def make_sing_b_dm(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
) -> NDArray[T]:
    r"""Make the single boson density matrix :math:`\langle b \rangle`.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.

    Returns:
        Single boson density matrix.
    """
    func, kwargs = self._load_function(
        "make_sing_b_dm",
        eris=eris,
        amplitudes=amplitudes,
        lambdas=lambdas,
    )
    res: NDArray[T] = func(**kwargs)
    return res

def make_rdm1_b(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
    unshifted: bool = True,
    hermitise: bool = True,
) -> NDArray[T]:
    r"""Make the one-particle boson reduced density matrix :math:`\langle b^+ c \rangle`.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.
        unshifted: If `self.options.shift` is `True`, return the unshifted density matrix. Has
            no effect if `self.options.shift` is `False`.
        hermitise: Hermitise the density matrix.

    Returns:
        One-particle boson reduced density matrix.
    """
    func, kwargs = self._load_function(
        "make_rdm1_b",
        eris=eris,
        amplitudes=amplitudes,
        lambdas=lambdas,
    )
    dm: NDArray[T] = func(**kwargs)

    if hermitise:
        dm = 0.5 * (dm + dm.T)

    if unshifted and self.options.shift:
        xi = self.xi
        dm_b = util.einsum("ni->i", self.make_sing_b_dm())
        dm -= util.einsum("ij,i->ij", np.eye(dm.shape[0]), xi * dm_b - xi**2.0)

    return dm

@abstractmethod
def make_rdm1_f(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
    hermitise: bool = True,
) -> Any:
    r"""Make the one-particle fermionic reduced density matrix :math:`\langle i^+ j \rangle`.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.
        hermitise: Hermitise the density matrix.

    Returns:
        One-particle fermion reduced density matrix.
    """
    pass

@abstractmethod
def make_rdm2_f(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
    hermitise: bool = True,
) -> Any:
    r"""Make the two-particle fermionic reduced density matrix :math:`\langle i^+j^+lk \rangle`.

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.
        hermitise: Hermitise the density matrix.

    Returns:
        Two-particle fermion reduced density matrix.
    """
    pass

@abstractmethod
def make_eb_coup_rdm(
    self,
    eris: Optional[ERIsInputType] = None,
    amplitudes: Optional[Namespace[SpinArrayType]] = None,
    lambdas: Optional[Namespace[SpinArrayType]] = None,
    unshifted: bool = True,
    hermitise: bool = True,
) -> Any:
    r"""Make the electron-boson coupling reduced density matrix.

    .. math::
        \langle b^+ i^+ j \rangle

    and

    .. math::
        \langle b i^+ j \rangle

    Args:
        eris: Electron repulsion integrals.
        amplitudes: Cluster amplitudes.
        lambdas: Cluster lambda amplitudes.
        unshifted: If `self.options.shift` is `True`, return the unshifted density matrix. Has
            no effect if `self.options.shift` is `False`.
        hermitise: Hermitise the density matrix.

    Returns:
        Electron-boson coupling reduced density matrix.
    """
    pass

@abstractmethod
def energy_sum(self, *args: str, signs_dict: Optional[dict[str, str]] = None) -> NDArray[T]:
    """Get a direct sum of energies.

    Args:
        *args: Energies to sum. Subclass should specify a subscript, and optionally spins.
        signs_dict: Signs of the energies in the sum. Default sets `("o", "O", "i")` to be
            positive, and `("v", "V", "a", "b")` to be negative.

    Returns:
        Sum of energies.
    """
    pass

@abstractmethod
def amplitudes_to_vector(self, amplitudes: Namespace[SpinArrayType]) -> NDArray[T]:
    """Construct a vector containing all of the amplitudes used in the given ansatz.

    Args:
        amplitudes: Cluster amplitudes.

    Returns:
        Cluster amplitudes as a vector.
    """
    pass

@abstractmethod
def vector_to_amplitudes(self, vector: NDArray[T]) -> Namespace[SpinArrayType]:
    """Construct a namespace of amplitudes from a vector.

    Args:
        vector: Cluster amplitudes as a vector.

    Returns:
        Cluster amplitudes.
    """
    pass

@abstractmethod
def lambdas_to_vector(self, lambdas: Namespace[SpinArrayType]) -> NDArray[T]:
    """Construct a vector containing all of the lambda amplitudes used in the given ansatz.

    Args:
        lambdas: Cluster lambda amplitudes.

    Returns:
        Cluster lambda amplitudes as a vector.
    """
    pass

@abstractmethod
def vector_to_lambdas(self, vector: NDArray[T]) -> Namespace[SpinArrayType]:
    """Construct a namespace of lambda amplitudes from a vector.

    Args:
        vector: Cluster lambda amplitudes as a vector.

    Returns:
        Cluster lambda amplitudes.
    """
    pass

@abstractmethod
def init_space(self) -> SpaceType:
    """Initialise the fermionic space.

    Returns:
        Fermionic space. All fermionic degrees of freedom are assumed to be correlated.
    """
    pass

@abstractmethod
def get_fock(self) -> BaseFock:
    """Get the Fock matrix.

    Returns:
        Fock matrix.
    """
    pass

@abstractmethod
def get_eris(self, eris: Optional[ERIsInputType] = None) -> BaseERIs:
    """Get the electron repulsion integrals.

    Args:
        eris: Input electron repulsion integrals.

    Returns:
        Electron repulsion integrals.
    """
    pass

def get_g(self) -> BaseElectronBoson:
    """Get the blocks of the electron-boson coupling matrix.

    This matrix corresponds to the bosonic annihilation operator.

    Returns:
        Electron-boson coupling matrix.
    """
    return self.ElectronBoson(self, array=self.bare_g)

@abstractmethod
def get_mean_field_G(self) -> Any:
    """Get the mean-field boson non-conserving term.

    Returns:
        Mean-field boson non-conserving term.
    """
    pass