diff --git a/ebcc/__init__.py b/ebcc/__init__.py
index 08859970..ce8e119d 100644
--- a/ebcc/__init__.py
+++ b/ebcc/__init__.py
@@ -36,89 +36,18 @@
 
 __version__ = "1.4.3"
 
-import logging
 import os
-import subprocess
 import sys
 
-# --- Logging:
-
-
-def output(self, msg, *args, **kwargs):
-    """Output a message at the `"OUTPUT"` level."""
-    if self.isEnabledFor(25):
-        self._log(25, msg, args, **kwargs)
-
-
-default_log = logging.getLogger(__name__)
-default_log.setLevel(logging.INFO)
-default_log.addHandler(logging.StreamHandler(sys.stderr))
-logging.addLevelName(25, "OUTPUT")
-logging.Logger.output = output
-
-
-class NullLogger(logging.Logger):
-    """A logger that does nothing."""
-
-    def __init__(self, *args, **kwargs):
-        super().__init__("null")
-
-    def _log(self, level, msg, args, **kwargs):
-        pass
-
-
-HEADER = """        _
-       | |
-   ___ | |__    ___   ___
-  / _ \| '_ \  / __| / __|
- |  __/| |_) || (__ | (__
-  \___||_.__/  \___| \___|
-%s"""
-
-
-def init_logging(log):
-    """Initialise the logging with a header."""
-
-    if globals().get("_EBCC_LOG_INITIALISED", False):
-        return
 
-    # Print header
-    header_size = max([len(line) for line in HEADER.split("\n")])
-    log.info(HEADER % (" " * (header_size - len(__version__)) + __version__))
-
-    # Print versions of dependencies and ebcc
-    def get_git_hash(directory):
-        git_directory = os.path.join(directory, ".git")
-        cmd = ["git", "--git-dir=%s" % git_directory, "rev-parse", "--short", "HEAD"]
-        try:
-            git_hash = subprocess.check_output(
-                cmd, universal_newlines=True, stderr=subprocess.STDOUT
-            ).rstrip()
-        except subprocess.CalledProcessError:
-            git_hash = "N/A"
-        return git_hash
-
-    import numpy
-    import pyscf
-
-    log.info("numpy:")
-    log.info(" > Version:  %s" % numpy.__version__)
-    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(numpy.__file__), "..")))
-
-    log.info("pyscf:")
-    log.info(" > Version:  %s" % pyscf.__version__)
-    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(pyscf.__file__), "..")))
+# --- Import NumPy here to allow drop-in replacements
 
-    log.info("ebcc:")
-    log.info(" > Version:  %s" % __version__)
-    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(__file__), "..")))
+import numpy
 
-    # Environment variables
-    log.info("OMP_NUM_THREADS = %s" % os.environ.get("OMP_NUM_THREADS", ""))
 
-    log.info("")
+# --- Logging:
 
-    globals()["_EBCC_LOG_INITIALISED"] = True
+from ebcc.logging import default_log, init_logging, NullLogger
 
 
 # --- Types of ansatz supporting by the EBCC solvers:
@@ -126,11 +55,6 @@ def get_git_hash(directory):
 METHOD_TYPES = ["MP", "CC", "LCC", "QCI", "QCC", "DC"]
 
 
-# --- Import NumPy here to allow drop-in replacements
-
-import numpy
-
-
 # --- General constructor:
 
 from ebcc.gebcc import GEBCC
@@ -189,6 +113,7 @@ def constructor(mf, *args, **kwargs):
 from ebcc.brueckner import BruecknerGEBCC, BruecknerREBCC, BruecknerUEBCC
 from ebcc.space import Space
 
+
 # --- List available methods:
 
 
diff --git a/ebcc/brueckner.py b/ebcc/brueckner.py
index ad6ee004..1bd1f0db 100644
--- a/ebcc/brueckner.py
+++ b/ebcc/brueckner.py
@@ -5,10 +5,11 @@
 import scipy.linalg
 from pyscf import lib
 
-from ebcc import NullLogger
+from ebcc import NullLogger, init_logging
 from ebcc import numpy as np
 from ebcc import util
 from ebcc.damping import DIIS
+from ebcc.logging import ANSI
 from ebcc.precision import types
 
 
@@ -78,12 +79,14 @@ def __init__(self, cc, log=None, options=None, **kwargs):
         self.converged = False
 
         # Logging:
-        cc.log.info("Brueckner options:")
-        cc.log.info(" > e_tol:  %s", options.e_tol)
-        cc.log.info(" > t_tol:  %s", options.t_tol)
-        cc.log.info(" > max_iter:  %s", options.max_iter)
-        cc.log.info(" > diis_space:  %s", options.diis_space)
-        cc.log.info(" > damping:  %s", options.damping)
+        init_logging(cc.log)
+        cc.log.info(f"\n{ANSI.B}{ANSI.U}{self.name}{ANSI.R}")
+        cc.log.info(f"{ANSI.B}Options{ANSI.R}:")
+        cc.log.info(f" > e_tol:  {ANSI.y}{self.options.e_tol}{ANSI.R}")
+        cc.log.info(f" > t_tol:  {ANSI.y}{self.options.t_tol}{ANSI.R}")
+        cc.log.info(f" > max_iter:  {ANSI.y}{self.options.max_iter}{ANSI.R}")
+        cc.log.info(f" > diis_space:  {ANSI.y}{self.options.diis_space}{ANSI.R}")
+        cc.log.info(f" > damping:  {ANSI.y}{self.options.damping}{ANSI.R}")
         cc.log.debug("")
 
     def get_rotation_matrix(self, u_tot=None, diis=None, t1=None):
@@ -292,17 +295,14 @@ def kernel(self):
         u_tot = None
 
         self.cc.log.output("Solving for Brueckner orbitals.")
-        self.cc.log.info(
-            "%4s %16s %18s %8s %13s %13s",
-            "Iter",
-            "Energy (corr.)",
-            "Energy (tot.)",
-            "Conv.",
-            "Δ(Energy)",
-            "|T1|",
+        self.cc.log.debug("")
+        self.log.info(
+            f"{ANSI.B}{'Iter':>4s} {'Energy (corr.)':>16s} {'Energy (tot.)':>18s} "
+            f"{'Conv.':>8s} {'Δ(Energy)':>13s} {'|T1|':>13s}{ANSI.R}"
         )
         self.log.info(
-            "%4d %16.10f %18.10f %8s", 0, self.cc.e_corr, self.cc.e_tot, self.cc.converged
+            f"{0:4d} {self.cc.e_corr:16.10f} {self.cc.e_tot:18.10f} "
+            f"{[ANSI.r, ANSI.g][self.cc.converged]}{self.cc.converged!r:>8}{ANSI.R}"
         )
 
         converged = False
@@ -336,23 +336,25 @@ def kernel(self):
             de = abs(e_prev - self.cc.e_tot)
             dt = self.get_t1_norm()
 
+            # Log the iteration:
+            converged_e = de < self.options.e_tol
+            converged_t = dt < self.options.t_tol
             self.log.info(
-                "%4d %16.10f %18.10f %8s %13.3e %13.3e",
-                niter,
-                self.cc.e_corr,
-                self.cc.e_tot,
-                self.cc.converged,
-                de,
-                dt,
+                f"{niter:4d} {self.cc.e_corr:16.10f} {self.cc.e_tot:18.10f}"
+                f" {[ANSI.r, ANSI.g][self.cc.converged]}{self.cc.converged!r:>8}{ANSI.R}"
+                f" {[ANSI.r, ANSI.g][converged_e]}{de:13.3e}{ANSI.R}"
+                f" {[ANSI.r, ANSI.g][converged_t]}{dt:13.3e}{ANSI.R}"
             )
 
             # Check for convergence:
-            converged = de < self.options.e_tol and dt < self.options.t_tol
+            converged = converged_e and converged_t
             if converged:
-                self.cc.log.output("Converged.")
+                self.log.debug("")
+                self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
                 break
         else:
-            self.cc.log.warning("Failed to converge.")
+            self.log.debug("")
+            self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}")
 
         self.cc.log.debug("")
         self.cc.log.output("E(corr) = %.10f", self.cc.e_corr)
@@ -360,10 +362,14 @@ def kernel(self):
         self.cc.log.debug("")
         self.cc.log.debug("Time elapsed: %s", timer.format_time(timer()))
         self.cc.log.debug("")
-        self.cc.log.debug("")
 
         return self.cc.e_corr
 
+    @property
+    def name(self):
+        """Get a string representation of the method name."""
+        return self.spin_type + "B" + self.cc.ansatz.name
+
     @property
     def spin_type(self):
         """Return the spin type."""
diff --git a/ebcc/logging.py b/ebcc/logging.py
new file mode 100644
index 00000000..b3561132
--- /dev/null
+++ b/ebcc/logging.py
@@ -0,0 +1,118 @@
+"""Logging."""
+
+import logging
+import os
+import subprocess
+import sys
+
+from ebcc import __version__
+from ebcc.util import Namespace
+
+HEADER = """        _
+       | |
+   ___ | |__    ___   ___
+  / _ \| '_ \  / __| / __|
+ |  __/| |_) || (__ | (__
+  \___||_.__/  \___| \___|
+%s"""  # noqa: W605
+
+
+def output(self, msg, *args, **kwargs):
+    """Output a message at the `"OUTPUT"` level."""
+    if self.isEnabledFor(25):
+        self._log(25, msg, args, **kwargs)
+
+
+default_log = logging.getLogger(__name__)
+default_log.setLevel(logging.INFO)
+default_log.addHandler(logging.StreamHandler(sys.stderr))
+logging.addLevelName(25, "OUTPUT")
+logging.Logger.output = output
+
+
+class NullLogger(logging.Logger):
+    """A logger that does nothing."""
+
+    def __init__(self, *args, **kwargs):
+        super().__init__("null")
+
+    def _log(self, level, msg, args, **kwargs):
+        pass
+
+
+def init_logging(log):
+    """Initialise the logging with a header."""
+
+    if globals().get("_EBCC_LOG_INITIALISED", False):
+        return
+
+    # Print header
+    header_size = max([len(line) for line in HEADER.split("\n")])
+    space = " " * (header_size - len(__version__))
+    log.info(f"{ANSI.B}{HEADER}{ANSI.R}" % f"{space}{ANSI.B}{__version__}{ANSI.R}")
+
+    # Print versions of dependencies and ebcc
+    def get_git_hash(directory):
+        git_directory = os.path.join(directory, ".git")
+        cmd = ["git", "--git-dir=%s" % git_directory, "rev-parse", "--short", "HEAD"]
+        try:
+            git_hash = subprocess.check_output(
+                cmd, universal_newlines=True, stderr=subprocess.STDOUT
+            ).rstrip()
+        except subprocess.CalledProcessError:
+            git_hash = "N/A"
+        return git_hash
+
+    import numpy
+    import pyscf
+
+    log.info("numpy:")
+    log.info(" > Version:  %s" % numpy.__version__)
+    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(numpy.__file__), "..")))
+
+    log.info("pyscf:")
+    log.info(" > Version:  %s" % pyscf.__version__)
+    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(pyscf.__file__), "..")))
+
+    log.info("ebcc:")
+    log.info(" > Version:  %s" % __version__)
+    log.info(" > Git hash: %s" % get_git_hash(os.path.join(os.path.dirname(__file__), "..")))
+
+    # Environment variables
+    log.info("OMP_NUM_THREADS = %s" % os.environ.get("OMP_NUM_THREADS", ""))
+
+    log.debug("")
+
+    globals()["_EBCC_LOG_INITIALISED"] = True
+
+
+def _check_output(*args, **kwargs):
+    """
+    Call a command. If the return code is non-zero, an empty `bytes`
+    object is returned.
+    """
+    try:
+        return subprocess.check_output(*args, **kwargs)
+    except subprocess.CalledProcessError:
+        return bytes()
+
+
+ANSI = Namespace(
+    B="\x1b[1m",
+    H="\x1b[3m",
+    R="\x1b[m\x0f",
+    U="\x1b[4m",
+    b="\x1b[34m",
+    c="\x1b[36m",
+    g="\x1b[32m",
+    k="\x1b[30m",
+    m="\x1b[35m",
+    r="\x1b[31m",
+    w="\x1b[37m",
+    y="\x1b[33m",
+)
+
+
+def colour(text, *cs):
+    """Colour a string."""
+    return f"{''.join([ANSI[c] for c in cs])}{text}{ANSI[None]}"
diff --git a/ebcc/rebcc.py b/ebcc/rebcc.py
index a69e0a88..dff7ee72 100644
--- a/ebcc/rebcc.py
+++ b/ebcc/rebcc.py
@@ -14,6 +14,7 @@
 from ebcc.dump import Dump
 from ebcc.eris import RERIs
 from ebcc.fock import RFock
+from ebcc.logging import ANSI
 from ebcc.precision import types
 from ebcc.space import Space
 
@@ -322,19 +323,19 @@ def __init__(
 
         # Logging:
         init_logging(self.log)
-        self.log.info("%s", self.name)
-        self.log.info("%s", "*" * len(self.name))
+        self.log.info(f"\n{ANSI.B}{ANSI.U}{self.name}{ANSI.R}")
+        self.log.debug(f"{ANSI.B}{'*' * len(self.name)}{ANSI.R}")
         self.log.debug("")
-        self.log.info("Options:")
-        self.log.info(" > e_tol:  %s", self.options.e_tol)
-        self.log.info(" > t_tol:  %s", self.options.t_tol)
-        self.log.info(" > max_iter:  %s", self.options.max_iter)
-        self.log.info(" > diis_space:  %s", self.options.diis_space)
-        self.log.info(" > damping:  %s", self.options.damping)
+        self.log.info(f"{ANSI.B}Options{ANSI.R}:")
+        self.log.info(f" > e_tol:  {ANSI.y}{self.options.e_tol}{ANSI.R}")
+        self.log.info(f" > t_tol:  {ANSI.y}{self.options.t_tol}{ANSI.R}")
+        self.log.info(f" > max_iter:  {ANSI.y}{self.options.max_iter}{ANSI.R}")
+        self.log.info(f" > diis_space:  {ANSI.y}{self.options.diis_space}{ANSI.R}")
+        self.log.info(f" > damping:  {ANSI.y}{self.options.damping}{ANSI.R}")
         self.log.debug("")
-        self.log.info("Ansatz: %s", self.ansatz)
+        self.log.info(f"{ANSI.B}Ansatz{ANSI.R}: {ANSI.m}{self.ansatz}{ANSI.R}")
         self.log.debug("")
-        self.log.info("Space: %s", self.space)
+        self.log.info(f"{ANSI.B}Space{ANSI.R}: {ANSI.m}{self.space}{ANSI.R}")
         self.log.debug("")
 
     def kernel(self, eris=None):
@@ -366,18 +367,15 @@ def kernel(self, eris=None):
             amplitudes = self.amplitudes
 
         # Get the initial energy:
-        e_cc = e_init = self.energy(amplitudes=amplitudes, eris=eris)
+        e_cc = self.energy(amplitudes=amplitudes, eris=eris)
 
         self.log.output("Solving for excitation amplitudes.")
+        self.log.debug("")
         self.log.info(
-            "%4s %16s %18s %13s %13s",
-            "Iter",
-            "Energy (corr.)",
-            "Energy (tot.)",
-            "Δ(Energy)",
-            "Δ(Ampl.)",
+            f"{ANSI.B}{'Iter':>4s} {'Energy (corr.)':>16s} {'Energy (tot.)':>18s} "
+            f"{'Δ(Energy)':>13s} {'Δ(Ampl.)':>13s}{ANSI.R}"
         )
-        self.log.info("%4d %16.10f %18.10f", 0, e_init, e_init + self.mf.e_tot)
+        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:
@@ -401,24 +399,32 @@ def kernel(self, eris=None):
                 e_cc = self.energy(amplitudes=amplitudes, eris=eris)
                 de = abs(e_prev - e_cc)
 
+                # Log the iteration:
+                converged_e = de < self.options.e_tol
+                converged_t = dt < self.options.t_tol
                 self.log.info(
-                    "%4d %16.10f %18.10f %13.3e %13.3e", niter, e_cc, e_cc + self.mf.e_tot, de, dt
+                    f"{niter:4d} {e_cc:16.10f} {e_cc + self.mf.e_tot:18.10f}"
+                    f" {[ANSI.r, ANSI.g][converged_e]}{de:13.3e}{ANSI.R}"
+                    f" {[ANSI.r, ANSI.g][converged_t]}{dt:13.3e}{ANSI.R}"
                 )
 
                 # Check for convergence:
-                converged = de < self.options.e_tol and dt < self.options.t_tol
+                converged = converged_e and converged_t
                 if converged:
-                    self.log.output("Converged.")
+                    self.log.debug("")
+                    self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
                     break
             else:
-                self.log.warning("Failed to converge.")
+                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("E(pert) = %.10f", e_pert)
+                self.log.info(f"E(pert) = {e_pert:.10f}")
 
         else:
             converged = True
@@ -429,12 +435,11 @@ def kernel(self, eris=None):
         self.converged = converged
 
         self.log.debug("")
-        self.log.output("E(corr) = %.10f", self.e_corr)
-        self.log.output("E(tot)  = %.10f", self.e_tot)
+        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("")
-        self.log.debug("")
 
         return e_cc
 
@@ -483,7 +488,8 @@ def solve_lambda(self, amplitudes=None, eris=None):
         diis.damping = self.options.damping
 
         self.log.output("Solving for de-excitation (lambda) amplitudes.")
-        self.log.info("%4s %13s", "Iter", "Δ(Ampl.)")
+        self.log.debug("")
+        self.log.info(f"{ANSI.B}{'Iter':>4s} {'Δ(Ampl.)':>13s}{ANSI.R}")
 
         converged = False
         for niter in range(1, self.options.max_iter + 1):
@@ -501,15 +507,18 @@ def solve_lambda(self, amplitudes=None, eris=None):
             lambdas = self.vector_to_lambdas(vector)
             dl = np.linalg.norm(vector - self.lambdas_to_vector(lambdas_prev), ord=np.inf)
 
-            self.log.info("%4d %13.3e", niter, dl)
+            # Log the iteration:
+            converged = dl < self.options.t_tol
+            self.log.info(f"{niter:4d} {[ANSI.r, ANSI.g][converged]}{dl:13.3e}{ANSI.R}")
 
             # Check for convergence:
-            converged = dl < self.options.t_tol
             if converged:
-                self.log.output("Converged.")
+                self.log.debug("")
+                self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
                 break
         else:
-            self.log.warning("Failed to converge.")
+            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()))
@@ -2050,16 +2059,16 @@ def boson_coupling_rank(self):
         """Get an integer representation of the boson coupling rank."""
         return self.ansatz.boson_coupling_rank
 
-    @property
-    def spin_type(self):
-        """Get a string represent of the spin channel."""
-        return "R"
-
     @property
     def name(self):
         """Get a string representation of the method name."""
         return self.spin_type + self.ansatz.name
 
+    @property
+    def spin_type(self):
+        """Get a string represent of the spin channel."""
+        return "R"
+
     @property
     def mo_coeff(self):
         """
diff --git a/ebcc/reom.py b/ebcc/reom.py
index 0d8c42f0..1c7ce752 100644
--- a/ebcc/reom.py
+++ b/ebcc/reom.py
@@ -7,6 +7,7 @@
 
 from ebcc import numpy as np
 from ebcc import util
+from ebcc.logging import ANSI
 from ebcc.precision import types
 
 
@@ -48,11 +49,7 @@ class REOM(EOM):
     Options = Options
 
     def __init__(self, ebcc, options=None, **kwargs):
-        self.ebcc = ebcc
-        self.space = ebcc.space
-        self.ansatz = ebcc.ansatz
-        self.log = ebcc.log
-
+        # Options:
         if options is None:
             options = self.Options()
         self.options = options
@@ -60,22 +57,30 @@ def __init__(self, ebcc, options=None, **kwargs):
             setattr(self.options, key, val)
         for key, val in self.options.__dict__.items():
             if val is util.Inherited:
-                setattr(self.options, key, getattr(self.ebcc.options, key))
+                setattr(self.options, key, getattr(ebcc.options, key))
 
-        self.log.info("%s", self.name)
-        self.log.info("%s", "*" * len(self.name))
-        self.log.debug("")
-        self.log.debug("Options:")
-        self.log.info(" > nroots:     %s", self.options.nroots)
-        self.log.info(" > e_tol:      %s", self.options.e_tol)
-        self.log.info(" > max_iter:   %s", self.options.max_iter)
-        self.log.info(" > max_space:  %s", self.options.max_space)
-        self.log.debug("")
+        # Parameters:
+        self.ebcc = ebcc
+        self.space = ebcc.space
+        self.ansatz = ebcc.ansatz
+        self.log = ebcc.log
 
+        # Attributes:
         self.converged = False
         self.e = None
         self.v = None
 
+        # Logging:
+        self.log.info(f"\n{ANSI.B}{ANSI.U}{self.name}{ANSI.R}")
+        self.log.debug(f"{ANSI.B}{'*' * len(self.name)}{ANSI.R}")
+        self.log.debug("")
+        self.log.info(f"{ANSI.B}Options{ANSI.R}:")
+        self.log.info(f" > nroots:  {ANSI.y}{self.options.nroots}{ANSI.R}")
+        self.log.info(f" > e_tol:  {ANSI.y}{self.options.e_tol}{ANSI.R}")
+        self.log.info(f" > max_iter:  {ANSI.y}{self.options.max_iter}{ANSI.R}")
+        self.log.info(f" > max_space:  {ANSI.y}{self.options.max_space}{ANSI.R}")
+        self.log.debug("")
+
     def amplitudes_to_vector(self, *amplitudes):
         """Convert the amplitudes to a vector."""
         raise NotImplementedError
@@ -163,23 +168,45 @@ def callback(self, envs):
     def _quasiparticle_weight(self, r1):
         return np.linalg.norm(r1) ** 2
 
-    def davidson(self, guesses=None):
-        """Solve the EOM Hamiltonian using the Davidson solver."""
+    def davidson(self, eris=None, guesses=None):
+        """Solve the EOM Hamiltonian using the Davidson solver.
+
+        Parameters
+        ----------
+        eris : ERIs, optional
+            Electronic repulsion integrals. Default value is generated
+            using `self.ebcc.get_eris()`.
+        guesses : list of np.ndarray, optional
+            Initial guesses for the roots. Default value is generated
+            using `self.get_guesses()`.
+
+        Returns
+        -------
+        e : np.ndarray
+            The energies of the roots.
+        """
+
+        # Start a timer:
+        timer = util.Timer()
+
+        # Get the ERIs:
+        eris = self.ebcc.get_eris(eris)
 
         self.log.output(
             "Solving for %s excitations using the Davidson solver.", self.excitation_type.upper()
         )
 
-        eris = self.ebcc.get_eris()
+        # Get the matrix-vector products and the diagonal:
         matvecs = lambda vs: [self.matvec(v, eris=eris) for v in vs]
         diag = self.diag(eris=eris)
 
+        # Get the guesses:
         if guesses is None:
             guesses = self.get_guesses(diag=diag)
 
+        # Solve the EOM Hamiltonian:
         nroots = min(len(guesses), self.options.nroots)
         pick = self.get_pick(guesses=guesses)
-
         converged, e, v = lib.davidson_nosym1(
             matvecs,
             guesses,
@@ -193,24 +220,35 @@ def davidson(self, guesses=None):
             verbose=0,
         )
 
+        # Check for convergence:
         if all(converged):
-            self.log.output("Converged.")
+            self.log.debug("")
+            self.log.output(f"{ANSI.g}Converged.{ANSI.R}")
         else:
-            self.log.warning("Failed to converge %d roots." % sum(not c for c in converged))
+            self.log.debug("")
+            self.log.warning(
+                f"{ANSI.r}Failed to converge {sum(not c for c in converged)} roots.{ANSI.R}"
+            )
 
-        self.log.debug("")
+        # Update attributes:
+        self.converged = converged
+        self.e = e
+        self.v = v
 
-        self.log.output("%4s %16s %16s", "Root", "Energy", "QP Weight")
-        for n, (en, vn) in enumerate(zip(e, v)):
+        self.log.debug("")
+        self.log.output(
+            f"{ANSI.B}{'Root':>4s} {'Energy':>16s} {'Weight':>13s} {'Conv.':>8s}{ANSI.R}"
+        )
+        for n, (en, vn, cn) in enumerate(zip(e, v, converged)):
             r1n = self.vector_to_amplitudes(vn)[0]
             qpwt = self._quasiparticle_weight(r1n)
-            self.log.output("%4d %16.10f %16.5g" % (n, en, qpwt))
+            self.log.output(
+                f"{n:>4d} {en:>16.10f} {qpwt:>13.5g} " f"{[ANSI.r, ANSI.g][cn]}{cn!r:>8s}{ANSI.R}"
+            )
 
         self.log.debug("")
-
-        self.converged = converged
-        self.e = e
-        self.v = v
+        self.log.debug("Time elapsed: %s", timer.format_time(timer()))
+        self.log.debug("")
 
         return e