diff --git a/.gitignore b/.gitignore
index c1b170cb2..cd73b1d7c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -128,3 +128,6 @@ Pipfile.lock
 
 #Containers
 *.sif
+
+#Visual Studio Code
+.vscode
diff --git a/doc/users/user_inp.rst b/doc/users/user_inp.rst
index 190e0e662..7893191bf 100644
--- a/doc/users/user_inp.rst
+++ b/doc/users/user_inp.rst
@@ -100,8 +100,10 @@ Available print levels are:
 - ``print_level=0`` prints mainly properties
 - ``print_level=1`` adds timings for individual steps
 - ``print_level=2`` adds memory and timing information on ``OrbitalVector`` level
-- ``print_level=3`` adds memory and timing information on ``Orbital`` level
-- ``print_level>10`` adds a *lot* more output from deep within MRCPP
+- ``print_level=3`` adds details for individual terms of the Fock operator
+- ``print_level=4`` adds memory and timing information on ``Orbital`` level
+- ``print_level>=5`` adds debug information at MRChem level
+- ``print_level>=10`` adds debug information at MRCPP level
 
 
 MPI
@@ -572,7 +574,7 @@ affect the convergence rate.
 
 For response calculations, the important convergence threshold is that of the
 orbitals, and by default this is set one order of magnitude higher than the
-overall ``world_prec``. 
+overall ``world_prec``.
 
 .. note::
 
diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py
index 8eadda3cb..c0bc4f1f0 100644
--- a/python/mrchem/helpers.py
+++ b/python/mrchem/helpers.py
@@ -195,10 +195,39 @@ def write_scf_solver(user_dict, wf_dict):
     if start_prec < 0.0:
         start_prec = final_prec
 
+    # TODO: Improve label determination. Make it easier to extent for future env additions
+    # TODO: Move determination to its own function to clean up this one
+    # TODO: Control precision of x,y,z strengths in label
+    # Determine the appropriate `Environment` printout in SCF 
+    # The relevant algorithms are implicit solvation and finite external fields
+    env = user_dict['Environment']
+    ext = user_dict['ExternalFields']
+
+    # Keep track of active environments
+    has_solvent = env['run_environment']
+    has_external_fields = len(ext['electric_field']) > 0
+
+    # If no external fields, then the list will be empty
+    # Need to catch the exception and store placeholders
+    try:
+        x, y, z = ext['electric_field']
+    except ValueError:
+        x, y, z = None, None, None  # Useless placeholders
+
+    # Labels to aggregate
+    label_sol = 'PCM'
+    label_ext = f'Electric field ({x}, {y}, {z})'
+
+    # Aggregate labels
+    labels = [label_sol, label_ext]
+    envs = [has_solvent, has_external_fields]
+    label_env = ' ; '.join([label for label, has_env in zip(labels, envs) if has_env])
+
     scf_dict = user_dict["SCF"]
     solver_dict = {
         "method": wf_dict["method_name"],
         "relativity": wf_dict["relativity_name"],
+        "environment": label_env,
         "kain": scf_dict["kain"],
         "max_iter": scf_dict["max_iter"],
         "rotation": scf_dict["rotation"],
@@ -211,6 +240,7 @@ def write_scf_solver(user_dict, wf_dict):
         "orbital_thrs": scf_dict["orbital_thrs"],
         "helmholtz_prec": user_dict["Precisions"]["helmholtz_prec"]
     }
+
     return solver_dict
 
 
@@ -446,5 +476,3 @@ def parse_wf_method(user_dict):
         "dft_funcs": dft_funcs
     }
     return wf_dict
-
-
diff --git a/src/chemistry/Molecule.cpp b/src/chemistry/Molecule.cpp
index d51db7374..8473c9691 100644
--- a/src/chemistry/Molecule.cpp
+++ b/src/chemistry/Molecule.cpp
@@ -208,6 +208,54 @@ void Molecule::printGeometry() const {
     mrcpp::print::separator(0, '=', 2);
 }
 
+/** @brief Pretty output of solvation cavity spheres */
+void Molecule::printCavity() {
+    // Collect relevant quantities
+    Cavity &cavity = getCavity();
+    std::vector<mrcpp::Coord<3>> coords = cavity.getCoordinates();
+    std::vector<double> radii = cavity.getRadii();
+
+    // Set widths
+    auto w0 = Printer::getWidth() - 1;
+    auto w1 = 5;
+    auto w2 = 12;
+    auto w3 = 2 * w0 / 9;
+    auto w4 = w0 - w1 - w2 - 3 * w3;
+
+    // Build table column headers
+    std::stringstream o_head;
+    o_head << std::setw(w1) << "N";
+    o_head << std::setw(w2) << "Radius";
+    o_head << std::string(w4 - 1, ' ') << ':';
+    o_head << std::setw(w3) << "x";
+    o_head << std::setw(w3) << "y";
+    o_head << std::setw(w3) << "z";
+
+    // Print
+    mrcpp::print::header(0, "Solvation Cavity");
+    print_utils::scalar(0, "Cavity width", cavity.getWidth(), "", 6);
+    mrcpp::print::separator(0, '-');
+    println(0, o_head.str());
+    mrcpp::print::separator(0, '-');
+    for (int i = 0; i < coords.size(); i++) {
+        mrcpp::Coord<3> coord = coords[i];
+        double r = radii[i];
+        double x = coord[0];
+        double y = coord[1];
+        double z = coord[2];
+
+        std::stringstream o_coord;
+        o_coord << std::setw(w1) << i;
+        o_coord << std::setw(w2) << std::setprecision(6) << std::fixed << r;
+        o_coord << std::string(w4 - 1, ' ') << ':';
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << x;
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << y;
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << z;
+        println(0, o_coord.str());
+    }
+    mrcpp::print::separator(0, '=', 2);
+}
+
 void Molecule::printEnergies(const std::string &txt) const {
     energy.print(txt);
     epsilon.print(txt);
diff --git a/src/chemistry/Molecule.h b/src/chemistry/Molecule.h
index 19e4727da..a778a5575 100644
--- a/src/chemistry/Molecule.h
+++ b/src/chemistry/Molecule.h
@@ -106,6 +106,7 @@ class Molecule final {
     void printGeometry() const;
     void printEnergies(const std::string &txt) const;
     void printProperties() const;
+    void printCavity();
 
     void initPerturbedOrbitals(bool dynamic);
     void initCavity(std::vector<mrcpp::Coord<3>> &coords, std::vector<double> &R, double slope);
diff --git a/src/chemistry/PhysicalConstants.cpp b/src/chemistry/PhysicalConstants.cpp
index 91ef82c7e..bf28b9299 100644
--- a/src/chemistry/PhysicalConstants.cpp
+++ b/src/chemistry/PhysicalConstants.cpp
@@ -58,9 +58,11 @@ PhysicalConstants &PhysicalConstants::Initialize(const json &constants) {
 
 /** @brief Pretty print physical constants */
 void PhysicalConstants::Print() {
-    mrcpp::print::header(0, "Physical Constants");
+    mrcpp::print::separator(0, '~');
+    print_utils::centered_text(0, "Physical Constants");
+    mrcpp::print::separator(0, '~');
     print_utils::json(0, constants_, true);
-    mrcpp::print::separator(0, '=', 2);
+    mrcpp::print::separator(0, '~', 2);
 }
 
 } // namespace mrchem
diff --git a/src/driver.cpp b/src/driver.cpp
index d766b677b..0fa626e81 100644
--- a/src/driver.cpp
+++ b/src/driver.cpp
@@ -141,6 +141,7 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) {
         auto xyz = coord["xyz"];
         nuclei.push_back(atom, xyz);
     }
+    mol.printGeometry();
 
     if (json_mol.contains("cavity")) {
         auto json_cavity = json_mol["cavity"];
@@ -153,8 +154,8 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) {
         auto width = json_cavity["width"];
 
         mol.initCavity(coords, radii, width);
+        // mol.printCavity();
     }
-    mol.printGeometry();
 }
 
 void driver::init_properties(const json &json_prop, Molecule &mol) {
@@ -222,7 +223,7 @@ void driver::init_properties(const json &json_prop, Molecule &mol) {
  * This function expects the "scf_calculation" subsection of the input.
  */
 json driver::scf::run(const json &json_scf, Molecule &mol) {
-    print_utils::headline(0, "Computing Ground State Wavefunction");
+    // print_utils::headline(0, "Computing Ground State Wavefunction");
     json json_out = {{"success", true}};
 
     if (json_scf.contains("properties")) driver::init_properties(json_scf["properties"], mol);
@@ -240,6 +241,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
     ///////////////////////////////////////////////////////////
     ///////////////   Setting Up Initial Guess   //////////////
     ///////////////////////////////////////////////////////////
+    print_utils::headline(0, "Computing Initial Guess Wavefunction");
     const auto &json_guess = json_scf["initial_guess"];
     if (scf::guess_orbitals(json_guess, mol)) {
         scf::guess_energy(json_guess, mol, F);
@@ -255,9 +257,12 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
 
     // Run GroundStateSolver if present in input JSON
     if (json_scf.contains("scf_solver")) {
+        print_utils::headline(0, "Computing Ground State Wavefunction");
+
         auto kain = json_scf["scf_solver"]["kain"];
         auto method = json_scf["scf_solver"]["method"];
         auto relativity = json_scf["scf_solver"]["relativity"];
+        auto environment = json_scf["scf_solver"]["environment"];
         auto max_iter = json_scf["scf_solver"]["max_iter"];
         auto rotation = json_scf["scf_solver"]["rotation"];
         auto localize = json_scf["scf_solver"]["localize"];
@@ -275,6 +280,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
         solver.setLocalize(localize);
         solver.setMethodName(method);
         solver.setRelativityName(relativity);
+        solver.setEnvironmentName(environment);
         solver.setCheckpoint(checkpoint);
         solver.setCheckpointFile(file_chk);
         solver.setMaxIterations(max_iter);
@@ -1032,9 +1038,10 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         auto formulation = json_fock["reaction_operator"]["formulation"];
         auto density_type = json_fock["reaction_operator"]["density_type"];
         Permittivity dielectric_func(*cavity_r, eps_in_r, eps_out_r, formulation);
+        dielectric_func.printParameters();
 
         auto scrf_p = std::make_unique<SCRF>(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type);
-        auto V_R = std::make_shared<ReactionOperator>(scrf_p, Phi_p);
+        auto V_R = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
         F.getReactionOperator() = V_R;
     }
     ///////////////////////////////////////////////////////////
diff --git a/src/environment/Cavity.h b/src/environment/Cavity.h
index e55888ee5..165febb62 100644
--- a/src/environment/Cavity.h
+++ b/src/environment/Cavity.h
@@ -50,6 +50,7 @@ class Cavity final : public mrcpp::RepresentableFunction<3> {
     auto getGradVector() const { return this->gradvector; }
     std::vector<mrcpp::Coord<3>> getCoordinates() const { return centers; } //!< Returns #centers.
     std::vector<double> getRadii() const { return radii; }                  //!< Returns #radii.
+    double getWidth() { return width; }                                     //!< Returns #width.
     friend class Permittivity;
 
 protected:
diff --git a/src/environment/Permittivity.cpp b/src/environment/Permittivity.cpp
index 33cf9e38d..e22ae6f0a 100644
--- a/src/environment/Permittivity.cpp
+++ b/src/environment/Permittivity.cpp
@@ -44,4 +44,54 @@ double Permittivity::evalf(const mrcpp::Coord<3> &r) const {
     }
 }
 
+void Permittivity::printParameters() {
+    // Collect relevant quantities
+    Cavity cavity = getCavity();
+    std::vector<mrcpp::Coord<3>> coords = cavity.getCoordinates();
+    std::vector<double> radii = cavity.getRadii();
+
+    // Set widths
+    auto w0 = mrcpp::Printer::getWidth() - 1;
+    auto w1 = 5;
+    auto w2 = 12;
+    auto w3 = 2 * w0 / 9;
+    auto w4 = w0 - w1 - w2 - 3 * w3;
+
+    // Build table column headers
+    std::stringstream o_head;
+    o_head << std::setw(w1) << "N";
+    o_head << std::setw(w2) << "Radius";
+    o_head << std::string(w4 - 1, ' ') << ':';
+    o_head << std::setw(w3) << "x";
+    o_head << std::setw(w3) << "y";
+    o_head << std::setw(w3) << "z";
+
+    // Print
+    mrcpp::print::header(0, "Solvation Cavity");
+    print_utils::scalar(0, "Cavity width", cavity.getWidth(), "", 6);
+    print_utils::scalar(0, "Dielec. Const. (in)", getEpsIn(), "", 6);
+    print_utils::scalar(0, "Dielec. Const. (out)", getEpsOut(), "", 6);
+    print_utils::text(0, "Formulation", getFormulation(), true);
+    mrcpp::print::separator(0, '-');
+    println(0, o_head.str());
+    mrcpp::print::separator(0, '-');
+    for (int i = 0; i < coords.size(); i++) {
+        mrcpp::Coord<3> coord = coords[i];
+        double r = radii[i];
+        double x = coord[0];
+        double y = coord[1];
+        double z = coord[2];
+
+        std::stringstream o_coord;
+        o_coord << std::setw(w1) << i;
+        o_coord << std::setw(w2) << std::setprecision(6) << std::fixed << r;
+        o_coord << std::string(w4 - 1, ' ') << ':';
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << x;
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << y;
+        o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << z;
+        println(0, o_coord.str());
+    }
+    mrcpp::print::separator(0, '=', 2);
+}
+
 } // namespace mrchem
diff --git a/src/environment/Permittivity.h b/src/environment/Permittivity.h
index a7adac7db..b2f72d5b3 100644
--- a/src/environment/Permittivity.h
+++ b/src/environment/Permittivity.h
@@ -26,7 +26,9 @@
 #pragma once
 
 #include "Cavity.h"
+#include "utils/print_utils.h"
 #include <MRCPP/MWFunctions>
+#include <MRCPP/Printer>
 
 namespace mrchem {
 /** @class Permittivity
@@ -83,6 +85,15 @@ class Permittivity final : public mrcpp::RepresentableFunction<3> {
     /** @brief Returns the value of #epsilon_out. */
     auto getEpsOut() const { return this->epsilon_out; }
 
+    /** @brief Returns the cavity */
+    Cavity getCavity() const { return this->cavity; }
+
+    /** @brief Returns the formulation */
+    std::string getFormulation() const { return this->formulation; }
+
+    /** @brief Print parameters */
+    void printParameters();
+
 private:
     bool inverse = false;    //!< State of #evalf
     double epsilon_in;       //!< Dielectric constant describing the permittivity of free space.
diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp
index fd1fab2d8..81eec02d8 100644
--- a/src/environment/SCRF.cpp
+++ b/src/environment/SCRF.cpp
@@ -26,6 +26,8 @@
 #include "SCRF.h"
 
 #include <MRCPP/MWOperators>
+#include <MRCPP/Printer>
+#include <MRCPP/Timer>
 
 #include "chemistry/PhysicalConstants.h"
 #include "chemistry/chemistry_utils.h"
@@ -35,6 +37,8 @@
 #include "scf_solver/KAIN.h"
 #include "utils/print_utils.h"
 
+#include <nlohmann/json.hpp>
+
 using mrcpp::Printer;
 using mrcpp::Timer;
 
@@ -62,6 +66,7 @@ SCRF::SCRF(Permittivity e,
         , max_iter(max_iter)
         , history(kain_hist)
         , apply_prec(orb_prec)
+        , conv_thrs(1.0)
         , mo_residual(-1.0)
         , epsilon(e)
         , rho_nuc(false)
@@ -87,6 +92,14 @@ void SCRF::clear() {
     this->rho_tot.free(NUMBER::Real);
 }
 
+double SCRF::setConvergenceThreshold(double prec) {
+    // converge_thrs should be in the interval [prec, 1.0]
+    this->conv_thrs = prec;
+    bool dynamic = (this->convergence_criterion == "dynamic");
+    if (dynamic and this->mo_residual > 10 * prec) this->conv_thrs = std::min(1.0, this->mo_residual);
+    return this->conv_thrs;
+}
+
 void SCRF::setDCavity() {
     mrcpp::FunctionTree<3> *dx_cavity = new mrcpp::FunctionTree<3>(*MRA);
     mrcpp::FunctionTree<3> *dy_cavity = new mrcpp::FunctionTree<3>(*MRA);
@@ -98,6 +111,7 @@ void SCRF::setDCavity() {
 }
 
 void SCRF::computeDensities(OrbitalVector &Phi) {
+    Timer timer;
     resetQMFunction(this->rho_tot);
     Density rho_el(false);
     density::compute(this->apply_prec, rho_el, Phi, DensityType::Total);
@@ -109,6 +123,7 @@ void SCRF::computeDensities(OrbitalVector &Phi) {
     } else {
         qmfunction::add(this->rho_tot, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); // probably change this into a vector
     }
+    print_utils::qmfunction(3, "Vacuum density", this->rho_tot, timer);
 }
 
 void SCRF::computeGamma(QMFunction &potential, QMFunction &out_gamma) {
@@ -159,18 +174,18 @@ void SCRF::accelerateConvergence(QMFunction &dfunc, QMFunction &func, KAIN &kain
 }
 
 void SCRF::nestedSCRF(QMFunction V_vac) {
-    print_utils::headline(2, "Calculating Reaction Potential");
-
     KAIN kain(this->history);
+    kain.setLocalPrintLevel(10);
 
-    // converge_thrs should be in the interval [apply_prec, 1.0]
-    double converge_thrs = this->apply_prec;
-    if (this->convergence_criterion == "dynamic" and this->mo_residual > 10 * this->apply_prec) converge_thrs = std::min(1.0, this->mo_residual);
+    mrcpp::print::separator(3, '-');
 
-    double update = 10.0;
-    for (int iter = 1; update >= converge_thrs && iter <= max_iter; iter++) {
+    double update = 10.0, norm = 1.0;
+    int iter = 1;
+    while (update >= this->conv_thrs && iter <= max_iter) {
+        Timer t_iter;
         // solve the poisson equation
         QMFunction Vr_np1 = solvePoissonEquation(this->gamma_n);
+        norm = Vr_np1.norm();
 
         // use a convergence accelerator
         resetQMFunction(this->dVr_n);
@@ -201,11 +216,10 @@ void SCRF::nestedSCRF(QMFunction V_vac) {
         }
 
         updateCurrentGamma(gamma_np1);
-
-        print_utils::text(2, "Update          ", print_utils::dbl_to_str(update, 5, true));
-        print_utils::text(2, "Microiteration  ", std::to_string(iter));
+        printConvergenceRow(iter, norm, update, t_iter.elapsed());
+        iter++;
     }
-    println(2, " Converged Reaction Potential!");
+    mrcpp::print::separator(3, '-');
     this->dVr_n.real().clear();
     this->dVr_n.real().setZero();
     this->dgamma_n.real().clear();
@@ -213,15 +227,42 @@ void SCRF::nestedSCRF(QMFunction V_vac) {
     kain.clear();
 }
 
+void SCRF::printConvergenceRow(int i, double norm, double update, double time) const {
+    auto pprec = Printer::getPrecision();
+    auto w0 = Printer::getWidth() - 1;
+    auto w1 = 9;
+    auto w2 = 2 * w0 / 9;
+    auto w3 = w0 - w1 - 2 * w2;
+
+    std::string time_unit = " sec";
+    if (time < 0.01) {
+        time = time * 1000.0;
+        time_unit = "  ms";
+    } else if (time > 60.0) {
+        time = time / 60.0;
+        time_unit = " min";
+    }
+
+    std::stringstream o_txt;
+    o_txt << " Iter " << std::setw(w1 - 6) << i;
+    o_txt << " : " << std::setw(w3 - 3) << std::setprecision(2 * pprec) << std::scientific << norm;
+    o_txt << std::setw(w2) << std::setprecision(pprec) << std::scientific << update;
+    o_txt << std::setw(w2 - 4) << std::setprecision(2) << std::fixed << time << time_unit;
+    println(3, o_txt.str());
+}
+
 QMFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) {
     this->apply_prec = prec;
     computeDensities(*Phi);
+    Timer t_vac;
     QMFunction V_vac;
     V_vac.alloc(NUMBER::Real);
     mrcpp::apply(this->apply_prec, V_vac.real(), *poisson, this->rho_tot.real());
+    print_utils::qmfunction(3, "Vacuum potential", V_vac, t_vac);
 
     // set up the zero-th iteration potential and gamma, so the first iteration gamma and potentials can be made
 
+    Timer t_gamma;
     if ((not this->Vr_n.hasReal()) or (not this->gamma_n.hasReal())) {
         QMFunction gamma_0;
         QMFunction V_tot;
@@ -248,8 +289,11 @@ QMFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) {
         qmfunction::deep_copy(this->gamma_n, temp_gamma_n);
         temp_gamma_n.free(NUMBER::Real);
     }
+    print_utils::qmfunction(3, "Initial gamma", this->gamma_n, t_gamma);
 
+    Timer t_scrf;
     if (this->algorithm == "scrf") nestedSCRF(V_vac);
+    print_utils::qmfunction(3, "Reaction potential", this->Vr_n, t_scrf);
     return this->Vr_n;
 }
 
@@ -290,4 +334,34 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) {
     gamma_np1.free(NUMBER::Real);
 }
 
+void SCRF::printParameters() const {
+    bool dynamic = (this->convergence_criterion == "dynamic");
+
+    std::stringstream o_iter;
+    if (this->max_iter > 0) {
+        o_iter << this->max_iter;
+    } else {
+        o_iter << "Off";
+    }
+
+    std::stringstream o_kain;
+    if (this->history > 0) {
+        o_kain << this->history;
+    } else {
+        o_kain << "Off";
+    }
+
+    nlohmann::json data = {
+        {"Method                ", "SCRF"},
+        {"Max iterations        ", max_iter},
+        {"KAIN solver           ", o_kain.str()},
+        {"Density type          ", density_type},
+        {"Dynamic threshold     ", (dynamic) ? "On" : "Off"},
+    };
+
+    mrcpp::print::separator(3, '~');
+    print_utils::json(3, data, false);
+    mrcpp::print::separator(3, '~');
+}
+
 } // namespace mrchem
diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h
index 6f1d20a4d..8c006b396 100644
--- a/src/environment/SCRF.h
+++ b/src/environment/SCRF.h
@@ -51,6 +51,8 @@ class SCRF final {
     ~SCRF();
     void UpdateExternalDensity(Density new_density) { this->rho_ext = new_density; }
 
+    double setConvergenceThreshold(double prec);
+
     QMFunction &getCurrentReactionPotential() { return this->Vr_n; }
     QMFunction &getPreviousReactionPotential() { return this->Vr_nm1; }
     QMFunction &getCurrentDifferenceReactionPotential() { return this->dVr_n; }
@@ -59,6 +61,8 @@ class SCRF final {
     QMFunction &getPreviousGamma() { return this->gamma_nm1; }
     QMFunction &getCurrentDifferenceGamma() { return this->dgamma_n; }
 
+    Permittivity &getPermittivity() { return this->epsilon; }
+
     void updateMOResidual(double const err_t) { this->mo_residual = err_t; }
 
     friend class ReactionPotential;
@@ -75,6 +79,7 @@ class SCRF final {
     int max_iter;
     int history;
     double apply_prec;
+    double conv_thrs;
     double mo_residual;
 
     Permittivity epsilon;
@@ -114,5 +119,8 @@ class SCRF final {
     void resetQMFunction(QMFunction &function);
     void updateCurrentReactionPotential(QMFunction &Vr_np1);
     void updateCurrentGamma(QMFunction &gamma_np1);
+
+    void printParameters() const;
+    void printConvergenceRow(int i, double norm, double update, double time) const;
 };
 } // namespace mrchem
diff --git a/src/initial_guess/core.cpp b/src/initial_guess/core.cpp
index 80372b630..3b2a37463 100644
--- a/src/initial_guess/core.cpp
+++ b/src/initial_guess/core.cpp
@@ -112,7 +112,6 @@ bool initial_guess::core::setup(OrbitalVector &Phi, double prec, const Nuclei &n
     t_lap.start();
     OrbitalVector Psi;
     initial_guess::core::project_ao(Psi, prec, nucs, zeta);
-    ComplexMatrix S_m12 = orbital::calc_lowdin_matrix(Psi);
     if (plevel == 1) mrcpp::print::time(1, "Projecting Hydrogen AOs", t_lap);
 
     p.setup(prec);
@@ -134,7 +133,7 @@ bool initial_guess::core::setup(OrbitalVector &Phi, double prec, const Nuclei &n
     V.clear();
     p.clear();
 
-    mrcpp::print::footer(2, t_lap, 2);
+    mrcpp::print::footer(2, t_tot, 2);
     if (plevel == 1) mrcpp::print::footer(1, t_tot, 2);
     return true;
 }
@@ -258,10 +257,12 @@ void initial_guess::core::rotate_orbitals(OrbitalVector &Psi, double prec, Compl
 ComplexMatrix initial_guess::core::diagonalize(OrbitalVector &Phi, MomentumOperator &p, RankZeroOperator &V) {
     Timer t1;
     ComplexMatrix S_m12 = orbital::calc_lowdin_matrix(Phi);
+    mrcpp::print::separator(2, '-');
     ComplexMatrix t_tilde = qmoperator::calc_kinetic_matrix(p, Phi, Phi);
     ComplexMatrix v_tilde = V(Phi, Phi);
     ComplexMatrix f_tilde = t_tilde + v_tilde;
     ComplexMatrix f = S_m12.adjoint() * f_tilde * S_m12;
+    mrcpp::print::separator(2, '-');
     mrcpp::print::time(1, "Computing Fock matrix", t1);
 
     Timer t2;
diff --git a/src/initial_guess/sad.cpp b/src/initial_guess/sad.cpp
index ba6709e3b..9f41cf110 100644
--- a/src/initial_guess/sad.cpp
+++ b/src/initial_guess/sad.cpp
@@ -115,11 +115,11 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::core::project_ao(Psi, prec, nucs, zeta);
     if (plevel == 1) mrcpp::print::time(1, "Projecting Hydrogen AOs", t_lap);
 
-    mrcpp::print::header(2, "Building Fock operator");
+    if (plevel == 2) mrcpp::print::header(2, "Building Fock operator");
     t_lap.start();
     p.setup(prec);
     V.setup(prec);
-    mrcpp::print::footer(2, t_lap, 2);
+    if (plevel == 2) mrcpp::print::footer(2, t_lap, 2);
     if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_lap);
 
     // Compute Fock matrix
@@ -138,7 +138,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     p.clear();
     V.clear();
 
-    mrcpp::print::footer(2, t_lap, 2);
+    mrcpp::print::footer(2, t_tot, 2);
     if (plevel == 1) mrcpp::print::footer(1, t_tot, 2);
     return true;
 }
@@ -194,12 +194,11 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     OrbitalVector Psi;
     initial_guess::gto::project_ao(Psi, prec, nucs);
     if (plevel == 1) mrcpp::print::time(1, "Projecting GTO AOs", t_lap);
-
-    mrcpp::print::header(2, "Building Fock operator");
+    if (plevel == 2) mrcpp::print::header(2, "Building Fock operator");
     t_lap.start();
     p.setup(prec);
     V.setup(prec);
-    mrcpp::print::footer(2, t_lap, 2);
+    if (plevel == 2) mrcpp::print::footer(2, t_lap, 2);
     if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_lap);
 
     // Compute Fock matrix
@@ -218,7 +217,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     p.clear();
     V.clear();
 
-    mrcpp::print::footer(2, t_lap, 2);
+    mrcpp::print::footer(2, t_tot, 2);
     if (plevel == 1) mrcpp::print::footer(1, t_tot, 2);
     return true;
 }
@@ -257,11 +256,10 @@ void initial_guess::sad::project_atomic_densities(double prec, Density &rho_tot,
     rho_loc.alloc(NUMBER::Real);
     rho_loc.real().setZero();
 
-    auto tot_nuc = 0.0;
-    auto tot_rho = 0.0;
-
     Timer t_loc;
-    for (int k = 0; k < nucs.size(); k++) {
+    auto N_nucs = nucs.size();
+    DoubleVector charges = DoubleVector::Zero(2 * N_nucs);
+    for (int k = 0; k < N_nucs; k++) {
         if (mpi::orb_rank != k % mpi::orb_size) continue;
 
         const std::string &sym = nucs[k].getElement().getSymbol();
@@ -273,24 +271,28 @@ void initial_guess::sad::project_atomic_densities(double prec, Density &rho_tot,
         rho_loc.add(1.0, rho_k);
         rho_loc.crop(crop_prec);
 
-        auto nuc_charge = nucs[k].getCharge();
-        auto rho_charge = rho_k.integrate().real();
-        tot_nuc += nuc_charge;
-        tot_rho += rho_charge;
+        charges[k] = nucs[k].getCharge();
+        charges[N_nucs + k] = rho_k.integrate().real();
+    }
+    t_loc.stop();
+
+    Timer t_com;
+    mpi::allreduce_vector(charges, mpi::comm_orb);
+    density::allreduce_density(prec, rho_tot, rho_loc);
+    t_com.stop();
 
+    for (int k = 0; k < N_nucs; k++) {
         std::stringstream o_row;
         o_row << std::setw(w1) << k;
-        o_row << std::setw(w2) << sym;
+        o_row << std::setw(w2) << nucs[k].getElement().getSymbol();
         o_row << std::setw(w4) << " ";
-        o_row << std::setw(w3) << print_utils::dbl_to_str(nuc_charge, 2 * pprec, false);
-        o_row << std::setw(w3) << print_utils::dbl_to_str(rho_charge, 2 * pprec, false);
+        o_row << std::setw(w3) << print_utils::dbl_to_str(charges[k], 2 * pprec, false);
+        o_row << std::setw(w3) << print_utils::dbl_to_str(charges[N_nucs + k], 2 * pprec, false);
         println(2, o_row.str());
     }
-    t_loc.stop();
 
-    Timer t_com;
-    density::allreduce_density(prec, rho_tot, rho_loc);
-    t_com.stop();
+    auto tot_nuc = charges.segment(0, N_nucs).sum();
+    auto tot_rho = charges.segment(N_nucs, N_nucs).sum();
 
     std::stringstream o_row;
     o_row << " Total charge";
diff --git a/src/mrdft/MRDFT.cpp b/src/mrdft/MRDFT.cpp
index 4a443f12c..0b7d03851 100644
--- a/src/mrdft/MRDFT.cpp
+++ b/src/mrdft/MRDFT.cpp
@@ -60,12 +60,22 @@ namespace mrdft {
  * out_vec[2] = v_xc_b (XC beta potential)
  */
 mrcpp::FunctionTreeVector<3> MRDFT::evaluate(mrcpp::FunctionTreeVector<3> &inp) {
-    mrcpp::Timer timer;
+    mrcpp::Timer t_tot, t_pre;
     grid().unify(inp);
     functional().preprocess(inp);
     mrcpp::FunctionTreeVector<3> xcInpVec = functional().setupXCInput();
     mrcpp::FunctionTreeVector<3> ctrInpVec = functional().setupCtrInput();
 
+    int inpNodes = 0;
+    int inpSize = 0;
+    for (auto i = 0; i < xcInpVec.size(); i++) {
+        auto &f_i = mrcpp::get_func(xcInpVec, i);
+        inpNodes += f_i.getNNodes();
+        inpSize += f_i.getSizeNodes();
+    }
+    mrcpp::print::tree(3, "Preprocess input", inpNodes, inpSize, t_pre.elapsed());
+
+    mrcpp::Timer t_eval;
     int nCoefs = mrcpp::get_func(inp, 0).getEndFuncNode(0).getNCoefs();
     int nOutCtr = functional().getCtrOutputLength();
     int nFcs = functional().getXCOutputLength();
@@ -135,22 +145,34 @@ mrcpp::FunctionTreeVector<3> MRDFT::evaluate(mrcpp::FunctionTreeVector<3> &inp)
     */
 
     // Reconstruct contracted output functions
-    int totNodes = 0;
-    int totSize = 0;
+    int ctrNodes = 0;
+    int ctrSize = 0;
     for (auto i = 0; i < ctrOutVec.size(); i++) {
         auto &f_i = mrcpp::get_func(ctrOutVec, i);
         f_i.mwTransform(mrcpp::BottomUp);
         f_i.calcSquareNorm();
         // TODO? insert a crop
-        totNodes += f_i.getNNodes();
-        totSize += f_i.getSizeNodes();
+        ctrNodes += f_i.getNNodes();
+        ctrSize += f_i.getSizeNodes();
     }
+    mrcpp::print::tree(3, "Evaluate functional", ctrNodes, ctrSize, t_eval.elapsed());
 
+    mrcpp::Timer t_post;
     auto potOutVec = functional().postprocess(ctrOutVec);
     mrcpp::clear(ctrOutVec, true);
     functional().clear();
-    auto t = timer.elapsed();
-    mrcpp::print::tree(2, "XC evaluate", totNodes, totSize, t);
+
+    int outNodes = 0;
+    int outSize = 0;
+    for (auto i = 0; i < potOutVec.size(); i++) {
+        auto &f_i = mrcpp::get_func(potOutVec, i);
+        f_i.mwTransform(mrcpp::BottomUp);
+        f_i.calcSquareNorm();
+        // TODO? insert a crop
+        outNodes += f_i.getNNodes();
+        outSize += f_i.getSizeNodes();
+    }
+    mrcpp::print::tree(3, "Postprocess potential", outNodes, outSize, t_post.elapsed());
 
     return potOutVec;
 }
diff --git a/src/parallel.cpp b/src/parallel.cpp
index 223f1d74b..370b19521 100644
--- a/src/parallel.cpp
+++ b/src/parallel.cpp
@@ -196,7 +196,7 @@ void mpi::initialize() {
 void mpi::finalize() {
 #ifdef MRCHEM_HAS_MPI
     if (mpi::bank_size > 0 and mpi::grand_master()) {
-        println(3, " max data in bank " << dataBank.get_maxtotalsize() << " MB ");
+        println(4, " max data in bank " << dataBank.get_maxtotalsize() << " MB ");
         dataBank.close();
     }
     MPI_Barrier(MPI_COMM_WORLD); // to ensure everybody got here
diff --git a/src/properties/SCFEnergy.h b/src/properties/SCFEnergy.h
index c82dd6b78..2a9a79b4e 100644
--- a/src/properties/SCFEnergy.h
+++ b/src/properties/SCFEnergy.h
@@ -74,6 +74,9 @@ class SCFEnergy final {
         auto E_kJ = E_au * PhysicalConstants::get("hartree2kjmol");
         auto E_kcal = E_au * PhysicalConstants::get("hartree2kcalmol");
 
+        bool has_ext = (std::abs(E_eext) > mrcpp::MachineZero) || (std::abs(E_next) > mrcpp::MachineZero);
+        bool has_react = (std::abs(Er_el) > mrcpp::MachineZero) || (std::abs(Er_nuc) > mrcpp::MachineZero);
+
         auto pprec = 2 * mrcpp::Printer::getPrecision();
         mrcpp::print::header(0, "Molecular Energy (" + id + ")");
         print_utils::scalar(0, "Kinetic energy   ", E_kin,  "(au)", pprec, false);
@@ -81,13 +84,19 @@ class SCFEnergy final {
         print_utils::scalar(0, "Coulomb energy   ", E_ee,   "(au)", pprec, false);
         print_utils::scalar(0, "Exchange energy  ", E_x,    "(au)", pprec, false);
         print_utils::scalar(0, "X-C energy       ", E_xc,   "(au)", pprec, false);
-        print_utils::scalar(0, "Ext. field (el)  ", E_eext, "(au)", pprec, false);
-        mrcpp::print::separator(0, '-');
         print_utils::scalar(0, "N-N energy       ", E_nn,   "(au)", pprec, false);
-        print_utils::scalar(0, "Ext. field (nuc) ", E_next, "(au)", pprec, false);
-        print_utils::scalar(0, "Tot. Reac. Energy", Er_tot,  "(au)", pprec, false);
-        print_utils::scalar(0, "El. Reac. Energy ", Er_el,  "(au)", pprec, false);
-        print_utils::scalar(0, "Nuc. Reac. Energy", Er_nuc,  "(au)", pprec, false);
+        if (has_ext) {
+            mrcpp::print::separator(0, '-');
+            print_utils::scalar(0, "External field (el)  ", E_eext, "(au)", pprec, false);
+            print_utils::scalar(0, "External field (nuc) ", E_next, "(au)", pprec, false);
+            print_utils::scalar(0, "External field (tot) ", E_eext + E_next, "(au)", pprec, false);
+        }
+        if (has_react) {
+            mrcpp::print::separator(0, '-');
+            print_utils::scalar(0, "Reaction energy (el)  ", Er_el,  "(au)", pprec, false);
+            print_utils::scalar(0, "Reaction energy (nuc) ", Er_nuc,  "(au)", pprec, false);
+            print_utils::scalar(0, "Reaction energy (tot) ", Er_tot,  "(au)", pprec, false);
+        }
         mrcpp::print::separator(0, '-');
         print_utils::scalar(0, "Electronic energy", E_el,   "(au)", pprec, false);
         print_utils::scalar(0, "Nuclear energy   ", E_nuc,  "(au)", pprec, false);
diff --git a/src/qmfunctions/orbital_utils.cpp b/src/qmfunctions/orbital_utils.cpp
index 97eff6b2d..6be524a51 100644
--- a/src/qmfunctions/orbital_utils.cpp
+++ b/src/qmfunctions/orbital_utils.cpp
@@ -1593,18 +1593,22 @@ int orbital::get_electron_number(const OrbitalVector &Phi, int spin) {
     return nElectrons;
 }
 
-/** @brief Returns the total number of nodes in the vector */
-int orbital::get_n_nodes(const OrbitalVector &Phi) {
-    int nNodes = 0;
-    for (const auto &phi_i : Phi) nNodes += phi_i.getNNodes(NUMBER::Total);
-    return nNodes;
+/** @brief Returns the total number of nodes in the vector, toggle to get average. */
+int orbital::get_n_nodes(const OrbitalVector &Phi, bool avg) {
+    long long totNodes = 0;
+    for (const auto &phi_i : Phi) totNodes += phi_i.getNNodes(NUMBER::Total);
+    if (avg) totNodes /= Phi.size();
+    if (totNodes > INT_MAX) MSG_WARN("Integer overflow: " << totNodes);
+    return static_cast<int>(totNodes);
 }
 
-/** @brief Returns the size of the coefficients of all nodes in the vector in kBytes */
-int orbital::get_size_nodes(const OrbitalVector &Phi) {
-    int tot_size = 0;
-    for (const auto &phi_i : Phi) tot_size += phi_i.getSizeNodes(NUMBER::Total);
-    return tot_size;
+/** @brief Returns the size of the coefficients of all nodes in the vector in kBytes, toggle to get average.*/
+int orbital::get_size_nodes(const OrbitalVector &Phi, bool avg) {
+    long long totSize = 0;
+    for (const auto &phi_i : Phi) totSize += phi_i.getSizeNodes(NUMBER::Total);
+    if (avg) totSize /= Phi.size();
+    if (totSize > INT_MAX) MSG_WARN("Integer overflow: " << totSize);
+    return static_cast<int>(totSize);
 }
 
 /** @brief Returns a vector containing the orbital spins */
diff --git a/src/qmfunctions/orbital_utils.h b/src/qmfunctions/orbital_utils.h
index a10db9adf..9e0efe2d5 100644
--- a/src/qmfunctions/orbital_utils.h
+++ b/src/qmfunctions/orbital_utils.h
@@ -81,8 +81,8 @@ int size_beta(const OrbitalVector &Phi);
 int get_multiplicity(const OrbitalVector &Phi);
 int get_electron_number(const OrbitalVector &Phi, int spin = SPIN::Paired);
 int start_index(const OrbitalVector &Phi, int spin);
-int get_n_nodes(const OrbitalVector &Phi);
-int get_size_nodes(const OrbitalVector &Phi);
+int get_n_nodes(const OrbitalVector &Phi, bool avg = false);
+int get_size_nodes(const OrbitalVector &Phi, bool avg = false);
 bool orbital_vector_is_sane(const OrbitalVector &Phi);
 
 void set_spins(OrbitalVector &Phi, const IntVector &spins);
diff --git a/src/qmoperators/one_electron/CMakeLists.txt b/src/qmoperators/one_electron/CMakeLists.txt
index 2291d40b7..264297579 100644
--- a/src/qmoperators/one_electron/CMakeLists.txt
+++ b/src/qmoperators/one_electron/CMakeLists.txt
@@ -1,3 +1,4 @@
 target_sources(mrchem PRIVATE
     ${CMAKE_CURRENT_SOURCE_DIR}/NuclearOperator.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/ZoraOperator.cpp
     )
diff --git a/src/qmoperators/one_electron/ZoraOperator.cpp b/src/qmoperators/one_electron/ZoraOperator.cpp
new file mode 100644
index 000000000..a31b071e3
--- /dev/null
+++ b/src/qmoperators/one_electron/ZoraOperator.cpp
@@ -0,0 +1,69 @@
+/*
+ * MRChem, a numerical real-space code for molecular electronic structure
+ * calculations within the self-consistent field (SCF) approximations of quantum
+ * chemistry (Hartree-Fock and Density Functional Theory).
+ * Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+ *
+ * This file is part of MRChem.
+ *
+ * MRChem is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MRChem is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with MRChem.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRChem, see:
+ * <https://mrchem.readthedocs.io/>
+ */
+
+#include "ZoraOperator.h"
+
+#include <MRCPP/Printer>
+#include <MRCPP/Timer>
+
+#include "qmfunctions/qmfunction_utils.h"
+#include "qmoperators/QMPotential.h"
+#include "utils/print_utils.h"
+
+using mrcpp::Printer;
+using mrcpp::Timer;
+
+namespace mrchem {
+
+ZoraOperator::ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse) {
+    Timer timer;
+    double two_cc = 2.0 * c * c;
+
+    auto k = std::make_shared<QMPotential>(1);
+    qmfunction::deep_copy(*k, vz);
+
+    if (k->hasImag()) MSG_ERROR("Inverse of complex function");
+    if (k->hasReal()) {
+        mrcpp::refine_grid(k->real(), 1);
+        if (inverse) {
+            k->real().map([two_cc](double val) { return (two_cc - val) / two_cc; });
+        } else {
+            k->real().map([two_cc](double val) { return two_cc / (two_cc - val); });
+        }
+        k->real().crop(proj_prec);
+    }
+
+    RankZeroOperator &kappa = (*this);
+    kappa = k;
+    if (inverse) {
+        kappa.name() = "kappa_m1";
+    } else {
+        kappa.name() = "kappa";
+    }
+    auto plevel = Printer::getPrintLevel();
+    print_utils::qmfunction(2, "ZORA operator (" + kappa.name() + ")", *k, timer);
+}
+
+} // namespace mrchem
diff --git a/src/qmoperators/one_electron/ZoraOperator.h b/src/qmoperators/one_electron/ZoraOperator.h
index d9abdbb73..0130b0c2e 100644
--- a/src/qmoperators/one_electron/ZoraOperator.h
+++ b/src/qmoperators/one_electron/ZoraOperator.h
@@ -25,39 +25,15 @@
 
 #pragma once
 
-#include "qmfunctions/qmfunction_utils.h"
-#include "qmoperators/QMPotential.h"
 #include "tensor/RankZeroOperator.h"
 
 namespace mrchem {
 
+class QMPotential;
+
 class ZoraOperator final : public RankZeroOperator {
 public:
-    ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse = false) {
-        double two_cc = 2.0 * c * c;
-
-        auto k = std::make_shared<QMPotential>(1);
-        qmfunction::deep_copy(*k, vz);
-
-        if (k->hasImag()) MSG_ERROR("Inverse of complex function");
-        if (k->hasReal()) {
-            mrcpp::refine_grid(k->real(), 1);
-            if (inverse) {
-                k->real().map([two_cc](double val) { return (two_cc - val) / two_cc; });
-            } else {
-                k->real().map([two_cc](double val) { return two_cc / (two_cc - val); });
-            }
-            k->real().crop(proj_prec);
-        }
-
-        RankZeroOperator &kappa = (*this);
-        kappa = k;
-        if (inverse) {
-            kappa.name() = "kappa_m1";
-        } else {
-            kappa.name() = "kappa";
-        }
-    }
+    ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse = false);
 };
 
 } // namespace mrchem
diff --git a/src/qmoperators/two_electron/CoulombPotential.cpp b/src/qmoperators/two_electron/CoulombPotential.cpp
index 4799944c5..fb9f016db 100644
--- a/src/qmoperators/two_electron/CoulombPotential.cpp
+++ b/src/qmoperators/two_electron/CoulombPotential.cpp
@@ -74,6 +74,11 @@ CoulombPotential::CoulombPotential(PoissonOperator_p P, OrbitalVector_p Phi, boo
 void CoulombPotential::setup(double prec) {
     if (isSetup(prec)) return;
     setApplyPrec(prec);
+    Timer timer;
+    auto plevel = Printer::getPrintLevel();
+    mrcpp::print::header(3, "Building Coulomb operator");
+    mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
+    mrcpp::print::separator(3, '-');
     if (hasDensity()) {
         setupGlobalPotential(prec);
     } else if (mpi::numerically_exact) {
@@ -86,6 +91,8 @@ void CoulombPotential::setup(double prec) {
         QMFunction V = setupLocalPotential(0.1 * prec);
         allreducePotential(0.1 * prec, V);
     }
+    if (plevel == 2) print_utils::qmfunction(2, "Coulomb operator", *this, timer);
+    mrcpp::print::footer(3, timer, 2);
 }
 
 /** @brief clear operator after application
@@ -124,7 +131,7 @@ void CoulombPotential::setupGlobalPotential(double prec) {
     V.alloc(NUMBER::Real);
     if (need_to_apply) mrcpp::apply(abs_prec, V.real(), P, rho.real());
     mpi::share_function(V, 0, 22445, mpi::comm_share);
-    print_utils::qmfunction(2, "Coulomb potential", V, timer);
+    print_utils::qmfunction(3, "Compute global potential", V, timer);
 }
 
 /** @brief compute Coulomb potential
@@ -148,7 +155,7 @@ QMFunction CoulombPotential::setupLocalPotential(double prec) {
     QMFunction V(false);
     V.alloc(NUMBER::Real);
     mrcpp::apply(abs_prec, V.real(), P, rho.real());
-    print_utils::qmfunction(2, "Coulomb potential", V, timer);
+    print_utils::qmfunction(3, "Compute local potential", V, timer);
 
     return V;
 }
@@ -183,7 +190,7 @@ void CoulombPotential::allreducePotential(double prec, QMFunction &V_loc) {
         mrcpp::copy_grid(V_tot.real(), V_loc.real());
         mrcpp::copy_func(V_tot.real(), V_loc.real());
     }
-    print_utils::qmfunction(2, "Allreduce Coulomb", V_tot, t_com);
+    print_utils::qmfunction(3, "Allreduce potential", V_tot, t_com);
 }
 
 } // namespace mrchem
diff --git a/src/qmoperators/two_electron/CoulombPotentialD1.cpp b/src/qmoperators/two_electron/CoulombPotentialD1.cpp
index 9f2ce960d..4989adb33 100644
--- a/src/qmoperators/two_electron/CoulombPotentialD1.cpp
+++ b/src/qmoperators/two_electron/CoulombPotentialD1.cpp
@@ -53,7 +53,7 @@ void CoulombPotentialD1::setupGlobalDensity(double prec) {
 
     Timer timer;
     density::compute(prec, rho, Phi, DensityType::Total);
-    print_utils::qmfunction(2, "Coulomb density", rho, timer);
+    print_utils::qmfunction(3, "Compute global density", rho, timer);
 }
 
 /** @brief compute local electron density (only own MPI orbitals)
@@ -71,7 +71,7 @@ void CoulombPotentialD1::setupLocalDensity(double prec) {
 
     Timer timer;
     density::compute_local(prec, rho, Phi, DensityType::Total);
-    print_utils::qmfunction(2, "Coulomb density", rho, timer);
+    print_utils::qmfunction(3, "Compute local density", rho, timer);
 }
 
 } // namespace mrchem
diff --git a/src/qmoperators/two_electron/CoulombPotentialD2.cpp b/src/qmoperators/two_electron/CoulombPotentialD2.cpp
index added2e14..ade5356a3 100644
--- a/src/qmoperators/two_electron/CoulombPotentialD2.cpp
+++ b/src/qmoperators/two_electron/CoulombPotentialD2.cpp
@@ -57,7 +57,7 @@ void CoulombPotentialD2::setupGlobalDensity(double prec) {
 
     Timer timer;
     density::compute(prec, rho, Phi, X, Y, DensityType::Total);
-    print_utils::qmfunction(2, "Coulomb density", rho, timer);
+    print_utils::qmfunction(3, "Compute global density", rho, timer);
 }
 
 void CoulombPotentialD2::setupLocalDensity(double prec) {
@@ -73,7 +73,7 @@ void CoulombPotentialD2::setupLocalDensity(double prec) {
 
     Timer timer;
     density::compute_local(prec, rho, Phi, X, Y, DensityType::Total);
-    print_utils::qmfunction(2, "Coulomb density", rho, timer);
+    print_utils::qmfunction(3, "Compute local density", rho, timer);
 }
 
 } // namespace mrchem
diff --git a/src/qmoperators/two_electron/ExchangePotential.cpp b/src/qmoperators/two_electron/ExchangePotential.cpp
index d8747e2b7..1fbfcf7be 100644
--- a/src/qmoperators/two_electron/ExchangePotential.cpp
+++ b/src/qmoperators/two_electron/ExchangePotential.cpp
@@ -117,10 +117,23 @@ double ExchangePotential::getSpinFactor(Orbital phi_i, Orbital phi_j) const {
  * and optionally compute the internal constributions.
  */
 void ExchangePotential::setup(double prec) {
+    Timer timer;
+    auto plevel = Printer::getPrintLevel();
+    mrcpp::print::header(3, "Building Exchange operator");
+    mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
+    mrcpp::print::separator(3, '-');
     if (mpi::world_size > 1 and mpi::bank_size < 1) MSG_ABORT("MPI bank required!");
     setApplyPrec(prec);
     setupBank();
     if (this->pre_compute) setupInternal(prec);
+
+    if (plevel == 2) {
+        auto t = timer.elapsed();
+        auto avgNodes = orbital::get_n_nodes(this->exchange, true);
+        auto avgSize = orbital::get_size_nodes(this->exchange, true);
+        mrcpp::print::tree(2, "Exchange operator (avg.)", avgNodes, avgSize, t);
+    }
+    mrcpp::print::footer(3, timer, 2);
 }
 
 /** @brief Clears the Exchange Operator
@@ -215,7 +228,7 @@ void ExchangePotential::calcExchange_kij(double prec, Orbital phi_k, Orbital phi
     }
     timer_jji.stop();
 
-    println(4,
+    println(5,
             " time " << (int)((float)timer_tot.elapsed() * 1000) << " ms "
                      << " mult1:" << (int)((float)timer_ij.elapsed() * 1000) << " Pot:" << (int)((float)timer_p.elapsed() * 1000) << " mult2:" << (int)((float)timer_kij.elapsed() * 1000) << " "
                      << (int)((float)timer_jji.elapsed() * 1000) << " Nnodes: " << N_i << " " << N_j << " " << N_ij << " " << N_p << " " << N_kij << " " << N_jji << " norms " << norm_ij << " "
diff --git a/src/qmoperators/two_electron/ExchangePotentialD1.cpp b/src/qmoperators/two_electron/ExchangePotentialD1.cpp
index 8a910ea4a..5c27c7768 100644
--- a/src/qmoperators/two_electron/ExchangePotentialD1.cpp
+++ b/src/qmoperators/two_electron/ExchangePotentialD1.cpp
@@ -64,7 +64,7 @@ void ExchangePotentialD1::setupBank() {
         if (mpi::my_orb(Phi[i])) PhiBank.put_orb(i, Phi[i]);
     }
     mpi::barrier(mpi::comm_orb);
-    mrcpp::print::time(4, "Setting up exchange bank", timer);
+    mrcpp::print::time(3, "Setting up exchange bank", timer);
 }
 
 /** @brief Clears rbital bank accounts.
@@ -116,10 +116,8 @@ Orbital ExchangePotentialD1::apply(Orbital phi_p) {
             MSG_WARN("Not computing exchange contributions that are not mine");
             return out_p;
         }
-        println(4, "On-the-fly exchange");
         return calcExchange(phi_p);
     } else {
-        println(4, "Precomputed exchange");
         return this->exchange[i];
     }
 }
@@ -144,7 +142,7 @@ QMOperatorVector ExchangePotentialD1::apply(QMOperator_p &O) {
  * The exchange potential is (pre)computed among the orbitals that define the operator
  */
 void ExchangePotentialD1::setupInternal(double prec) {
-    Timer timerT, timerS(false), timerR(false), t_calc(false), t_add(false), t_get(false), t_wait(false);
+    Timer t_tot, t_snd(false), t_orb(false), t_calc(false), t_add(false), t_get(false), t_wait(false);
     setApplyPrec(prec);
     if (this->exchange.size() != 0) MSG_ERROR("Exchange not properly cleared");
 
@@ -157,7 +155,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
     prec = mpi::numerically_exact ? -1.0 : prec;
     precf /= std::sqrt(1 * Phi.size());
     // Initialize this->exchange and compute own diagonal elements
-    Timer timerD;
+    Timer t_diag;
     for (auto &phi_i : Phi) {
         Orbital ex_iii = phi_i.paramCopy();
         t_calc.resume();
@@ -165,8 +163,9 @@ void ExchangePotentialD1::setupInternal(double prec) {
         t_calc.stop();
         Ex.push_back(ex_iii);
     }
-    mrcpp::print::time(4, "Exchange time diagonal", timerD);
+    t_diag.stop();
 
+    Timer t_offd;
     // We divide all the exchange contributions into a fixed number of tasks.
     // all "j" orbitals are fetched and stored, and used together with one "i" orbital
     // At the end of a task, the Exchange for j is summed up with the contributions
@@ -250,25 +249,25 @@ void ExchangePotentialD1::setupInternal(double prec) {
             i0 = iorb;
             Orbital phi_i;
 
-            timerR.resume();
+            t_orb.resume();
             if (bank_size > 0) {
                 PhiBank.get_orb(iorb, phi_i, 1); // fetch also own orbitals (simpler for clean up, and they are few)
                 iorb_vec.push_back(phi_i);
             } else {
                 iorb_vec.push_back(Phi[iorb]);
             }
-            timerR.stop();
+            t_orb.stop();
         }
 
         for (int j = 0; j < jtasks[task].size(); j++) {
             int jorb = jtasks[task][j];
             Orbital phi_j;
-            timerR.resume();
+            t_orb.resume();
             if (bank_size > 0)
                 PhiBank.get_orb(jorb, phi_j, 1);
             else
                 phi_j = Phi[jorb];
-            timerR.stop();
+            t_orb.stop();
             QMFunctionVector iijfunc_vec;
             ComplexVector coef_vec(N);
             for (int i = 0; i < iorb_vec.size(); i++) {
@@ -284,7 +283,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
                 calcExchange_kij(precf, phi_i, phi_i, phi_j, ex_iij, &ex_jji);
                 t_calc.stop();
                 if (ex_iij.norm() > prec) coef_vec[iijfunc_vec.size()] = j_fac;
-                timerS.resume();
+                t_snd.resume();
                 if (bank_size > 0) {
                     // store ex_jji
                     if (ex_iij.norm() > prec) iijfunc_vec.push_back(ex_iij);
@@ -295,7 +294,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
                     Ex[jorb].add(j_fac, ex_iij);
                 }
                 ex_jji.free(NUMBER::Total);
-                timerS.stop();
+                t_snd.stop();
             }
             Timer timerx;
             // fetch ready contributions to ex_j from others
@@ -320,14 +319,14 @@ void ExchangePotentialD1::setupInternal(double prec) {
                 // ex_j is sent to Bank
                 if (ex_j.hasReal() or ex_j.hasImag()) {
                     auto tT = timerx.elapsed();
-                    timerS.resume();
+                    t_snd.resume();
                     ex_j.crop(prec);
                     if (ex_j.norm() > prec) {
                         if (i0 < 0) MSG_ERROR("no exchange contributions?");
                         ExBank.put_orb(jorb + (i0 + N) * N, ex_j);
                         tasksMaster.put_readytask(jorb, i0 + N);
                     }
-                    timerS.stop();
+                    t_snd.stop();
                     ex_j.free(NUMBER::Total);
                 } else if (iijfunc_vec.size() > 0) {
                     MSG_ERROR("Exchange exists but has no real and no Imag parts");
@@ -342,7 +341,6 @@ void ExchangePotentialD1::setupInternal(double prec) {
     mpi::barrier(mpi::comm_orb);
     t_wait.stop();
 
-    IntVector sizes = IntVector::Zero(2 * N);
     for (int j = 0; j < N; j++) {
         if (not mpi::my_orb(Phi[j]) or bank_size == 0) continue; // fetch only own j
         std::vector<int> iVec = tasksMaster.get_readytask(j, 1);
@@ -385,26 +383,24 @@ void ExchangePotentialD1::setupInternal(double prec) {
             Ex[j].crop(prec);
             t_add.stop();
         }
-        sizes[j] = Ex[j].getNNodes(NUMBER::Total);
-        sizes[j + N] = Ex[j].getSizeNodes(NUMBER::Total);
     }
-    mrcpp::print::time(4, "Time rcv orbitals", timerR);
-    mrcpp::print::time(4, "Time send exchanges", timerS);
-    mrcpp::print::time(4, "Time rcv exchanges", t_get);
-    mrcpp::print::time(4, "Time wait others finished", t_wait);
-    mrcpp::print::time(4, "Time add exchanges", t_add);
-    mrcpp::print::time(4, "Time calculate exchanges", t_calc);
-
-    auto t = timerT.elapsed();
-    mpi::allreduce_vector(sizes, mpi::comm_orb);
-    long long nsum = 0;
-    for (int j = 0; j < N; j++) nsum += sizes[j];
-    long long msum = 0;
-    for (int j = 0; j < N; j++) msum += sizes[j + N];
-    int n = nsum / N;
-    int m = msum / N;
-    mrcpp::print::tree(2, "HF exchange (av.)", n, m, t);
-} // namespace mrchem
+    t_offd.stop();
+    mrcpp::print::time(3, "Time receiving orbitals", t_orb);
+    mrcpp::print::time(3, "Time receiving exchanges", t_get);
+    mrcpp::print::time(3, "Time sending exchanges", t_snd);
+    mrcpp::print::time(3, "Time adding exchanges", t_add);
+    mrcpp::print::time(3, "Time waiting for others", t_wait);
+    mrcpp::print::time(3, "Time computing exchanges", t_calc);
+    mrcpp::print::separator(3, '-');
+    mrcpp::print::time(3, "Time diagonal terms", t_diag);
+    mrcpp::print::time(3, "Time off-diagonal terms", t_offd);
+    mrcpp::print::separator(3, '-');
+
+    auto t = t_tot.elapsed() / N;
+    auto n = orbital::get_n_nodes(this->exchange, true);
+    auto m = orbital::get_size_nodes(this->exchange, true);
+    mrcpp::print::tree(3, "Average exchange term", n, m, t);
+}
 
 /** @brief Computes the exchange potential on the fly
  *
@@ -444,7 +440,7 @@ Orbital ExchangePotentialD1::calcExchange(Orbital phi_p) {
     Orbital ex_p = phi_p.paramCopy();
     Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     qmfunction::linear_combination(ex_p, coefs, func_vec, prec);
-    print_utils::qmfunction(3, "Applied exchange", ex_p, timer);
+    print_utils::qmfunction(4, "Applied exchange", ex_p, timer);
     return ex_p;
 }
 
diff --git a/src/qmoperators/two_electron/ExchangePotentialD2.cpp b/src/qmoperators/two_electron/ExchangePotentialD2.cpp
index 4533ccbfa..169719967 100644
--- a/src/qmoperators/two_electron/ExchangePotentialD2.cpp
+++ b/src/qmoperators/two_electron/ExchangePotentialD2.cpp
@@ -76,7 +76,7 @@ void ExchangePotentialD2::setupBank() {
         if (mpi::my_orb(Y[i])) YBank.put_orb(i, Y[i]);
     }
     mpi::barrier(mpi::comm_orb);
-    mrcpp::print::time(4, "Setting up exchange bank", timer);
+    mrcpp::print::time(3, "Setting up exchange bank", timer);
 }
 
 /** @brief Clears the Exchange Operator
@@ -144,7 +144,7 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) {
     Orbital out_p = phi_p.paramCopy();
     Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     qmfunction::linear_combination(out_p, coefs, func_vec, prec);
-    print_utils::qmfunction(3, "Applied exchange", out_p, timer);
+    print_utils::qmfunction(4, "Applied exchange", out_p, timer);
     return out_p;
 }
 
@@ -203,7 +203,7 @@ Orbital ExchangePotentialD2::dagger(Orbital phi_p) {
     Orbital ex_p = phi_p.paramCopy();
     Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     qmfunction::linear_combination(ex_p, coefs, func_vec, prec);
-    print_utils::qmfunction(3, "Applied exchange", ex_p, timer);
+    print_utils::qmfunction(4, "Applied exchange", ex_p, timer);
     return ex_p;
 }
 
diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp
index 4e6b9cdd9..91a03a93b 100644
--- a/src/qmoperators/two_electron/FockBuilder.cpp
+++ b/src/qmoperators/two_electron/FockBuilder.cpp
@@ -75,15 +75,22 @@ void FockBuilder::build(double exx) {
 void FockBuilder::setup(double prec) {
     Timer t_tot;
     auto plevel = Printer::getPrintLevel();
-    mrcpp::print::header(2, "Building Fock operator");
-    mrcpp::print::value(2, "Precision", prec, "(rel)", 5);
-    mrcpp::print::separator(2, '-');
+    if (plevel == 2) {
+        mrcpp::print::header(2, "Building Fock operator");
+        mrcpp::print::value(2, "Precision", prec, "(rel)", 5);
+        mrcpp::print::separator(2, '-');
+    }
     if (this->mom != nullptr) this->momentum().setup(prec);
     this->potential().setup(prec);
     this->perturbation().setup(prec);
 
     if (isZora()) {
+        Timer t_zora;
         auto c = getLightSpeed();
+        mrcpp::print::header(3, "Building ZORA operators");
+        mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
+        mrcpp::print::value(3, "Light speed", c, "(au)", 5);
+        mrcpp::print::separator(3, '-');
         auto vz = collectZoraBasePotential();
         this->kappa = std::make_shared<ZoraOperator>(*vz, c, prec, false);
         this->kappa_inv = std::make_shared<ZoraOperator>(*vz, c, prec, true);
@@ -91,10 +98,11 @@ void FockBuilder::setup(double prec) {
         this->kappa->setup(prec);
         this->kappa_inv->setup(prec);
         this->zora_base.setup(prec);
+        mrcpp::print::footer(3, t_zora, 2);
     };
 
     t_tot.stop();
-    mrcpp::print::footer(2, t_tot, 2);
+    if (plevel == 2) mrcpp::print::footer(2, t_tot, 2);
     if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_tot);
 }
 
@@ -203,12 +211,13 @@ ComplexMatrix FockBuilder::operator()(OrbitalVector &bra, OrbitalVector &ket) {
 }
 
 OrbitalVector FockBuilder::buildHelmholtzArgument(double prec, OrbitalVector Phi, ComplexMatrix F_mat, ComplexMatrix L_mat) {
-    Timer t_arg;
+    Timer t_tot;
     auto plevel = Printer::getPrintLevel();
     mrcpp::print::header(2, "Computing Helmholtz argument");
 
+    Timer t_rot;
     OrbitalVector Psi = orbital::rotate(Phi, L_mat - F_mat, prec);
-    mrcpp::print::time(2, "Rotating orbitals", t_arg);
+    mrcpp::print::time(2, "Rotating orbitals", t_rot);
 
     OrbitalVector out;
     if (isZora()) {
@@ -218,8 +227,8 @@ OrbitalVector FockBuilder::buildHelmholtzArgument(double prec, OrbitalVector Phi
     }
     Psi.clear();
 
-    mrcpp::print::footer(2, t_arg, 2);
-    if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_arg);
+    mrcpp::print::footer(2, t_tot, 2);
+    if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_tot);
     return out;
 }
 
@@ -239,17 +248,24 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     operOne.setup(prec);
     operThree.setup(prec);
 
+    // Compute OrbitalVectors
+    Timer t_1;
+    OrbitalVector termOne = operOne(Phi);
+    mrcpp::print::time(2, "Computing gradient term", t_1);
+
+    Timer t_2;
+    OrbitalVector termTwo = V(Phi);
+    mrcpp::print::time(2, "Computing potential term", t_2);
+
     // Compute transformed orbitals scaled by diagonal Fock elements
+    Timer t_3;
     OrbitalVector epsPhi = orbital::deep_copy(Phi);
     for (int i = 0; i < epsPhi.size(); i++) {
         if (not mpi::my_orb(epsPhi[i])) continue;
         epsPhi[i].rescale(eps[i] / two_cc);
     }
-
-    // Compute OrbitalVectors
-    OrbitalVector termOne = operOne(Phi);
-    OrbitalVector termTwo = V(Phi);
     OrbitalVector termThree = operThree(epsPhi);
+    mrcpp::print::time(2, "Computing rescaled potential term", t_3);
 
     auto normsOne = orbital::get_norms(termOne);
     auto normsTwo = orbital::get_norms(termTwo);
@@ -257,17 +273,23 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     auto normsPsi = orbital::get_norms(Psi);
 
     // Add up all the terms
-    OrbitalVector out = orbital::deep_copy(termOne);
-    for (int i = 0; i < out.size(); i++) {
-        if (not mpi::my_orb(out[i])) continue;
-        out[i].add(1.0, termTwo[i]);
-        out[i].add(1.0, termThree[i]);
-        out[i].add(1.0, Psi[i]);
+    Timer t_add;
+    OrbitalVector arg = orbital::deep_copy(termOne);
+    for (int i = 0; i < arg.size(); i++) {
+        if (not mpi::my_orb(arg[i])) continue;
+        arg[i].add(1.0, termTwo[i]);
+        arg[i].add(1.0, termThree[i]);
+        arg[i].add(1.0, Psi[i]);
     }
+    mrcpp::print::time(2, "Adding contributions", t_add);
 
     operThree.clear();
     operOne.clear();
-    return kappa_m1(out);
+
+    Timer t_kappa;
+    auto out = kappa_m1(arg);
+    mrcpp::print::time(2, "Applying kappa inverse", t_kappa);
+    return out;
 }
 
 // Non-relativistic Helmholtz argument
@@ -276,14 +298,18 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentNREL(OrbitalVector &Phi, Orbita
     RankZeroOperator &V = this->potential();
 
     // Compute OrbitalVectors
+    Timer t_pot;
     OrbitalVector termOne = V(Phi);
+    mrcpp::print::time(2, "Computing potential term", t_pot);
 
     // Add up all the terms
+    Timer t_add;
     OrbitalVector out = orbital::deep_copy(termOne);
     for (int i = 0; i < out.size(); i++) {
         if (not mpi::my_orb(out[i])) continue;
         out[i].add(1.0, Psi[i]);
     };
+    mrcpp::print::time(2, "Adding contributions", t_add);
     return out;
 }
 
@@ -294,6 +320,7 @@ void FockBuilder::setZoraType(bool has_nuc, bool has_coul, bool has_xc) {
 }
 
 std::shared_ptr<QMPotential> FockBuilder::collectZoraBasePotential() {
+    Timer timer;
     auto vz = std::make_shared<QMPotential>(1, false);
     if (zora_has_nuc) {
         if (getNuclearOperator() != nullptr) {
@@ -324,6 +351,7 @@ std::shared_ptr<QMPotential> FockBuilder::collectZoraBasePotential() {
             MSG_ERROR("ZORA: XC requested but not available");
         }
     }
+    print_utils::qmfunction(2, "ZORA operator (base)", *vz, timer);
     return vz;
 }
 
diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h
index 58f12a2e7..dae73f5f2 100644
--- a/src/qmoperators/two_electron/ReactionOperator.h
+++ b/src/qmoperators/two_electron/ReactionOperator.h
@@ -42,11 +42,12 @@ class SCRF;
 
 class ReactionOperator final : public RankZeroOperator {
 public:
-    ReactionOperator(std::unique_ptr<SCRF> &scrf_p, std::shared_ptr<mrchem::OrbitalVector> Phi_p) {
-        potential = std::make_shared<ReactionPotential>(scrf_p, Phi_p);
+    ReactionOperator(std::unique_ptr<SCRF> scrf_p, std::shared_ptr<mrchem::OrbitalVector> Phi_p) {
+        potential = std::make_shared<ReactionPotential>(std::move(scrf_p), Phi_p);
         // Invoke operator= to assign *this operator
         RankZeroOperator &V = (*this);
         V = potential;
+        V.name() = "V_r";
     }
 
     ComplexDouble trace(OrbitalVector &Phi) { return RankZeroOperator::trace(Phi); }
diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp
index da723d7b3..15f250693 100644
--- a/src/qmoperators/two_electron/ReactionPotential.cpp
+++ b/src/qmoperators/two_electron/ReactionPotential.cpp
@@ -23,6 +23,9 @@
  * <https://mrchem.readthedocs.io/>
  */
 
+#include <MRCPP/Printer>
+#include <MRCPP/Timer>
+
 #include "ReactionPotential.h"
 
 #include "qmfunctions/qmfunction_utils.h"
@@ -32,7 +35,7 @@ using OrbitalVector_p = std::shared_ptr<mrchem::OrbitalVector>;
 
 namespace mrchem {
 
-ReactionPotential::ReactionPotential(SCRF_p &scrf_p, OrbitalVector_p Phi_p)
+ReactionPotential::ReactionPotential(SCRF_p scrf_p, OrbitalVector_p Phi_p)
         : QMPotential(1, false)
         , helper(std::move(scrf_p))
         , Phi(Phi_p) {}
@@ -43,8 +46,19 @@ void ReactionPotential::setup(double prec) {
         this->first_iteration = false;
         return;
     }
+    auto thrs = this->helper->setConvergenceThreshold(prec);
+    mrcpp::Timer timer;
+    auto plevel = mrcpp::Printer::getPrintLevel();
+    mrcpp::print::separator(3, '=');
+    print_utils::centered_text(3, "Building Reaction operator");
+    this->helper->printParameters();
+    mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
+    mrcpp::print::value(3, "Threshold", thrs, "(abs)", 5);
+    mrcpp::print::separator(3, '-');
     auto potential = this->helper->setup(prec, this->Phi);
     qmfunction::deep_copy(*this, potential);
+    if (plevel == 2) print_utils::qmfunction(2, "Reaction operator", *this, timer);
+    mrcpp::print::footer(3, timer, 2);
 }
 
 void ReactionPotential::clear() {
diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h
index 4af12b514..7aade0014 100644
--- a/src/qmoperators/two_electron/ReactionPotential.h
+++ b/src/qmoperators/two_electron/ReactionPotential.h
@@ -44,7 +44,7 @@ class ReactionPotential final : public QMPotential {
     /** @brief Initializes the ReactionPotential class.
      *  @param scrf_p A SCRF instance which contains the parameters needed to compute the ReactionPotential.
      *  @param Phi_p A pointer to a vector which contains the orbitals optimized in the SCF procedure. */
-    ReactionPotential(std::unique_ptr<SCRF> &scrf_p, std::shared_ptr<mrchem::OrbitalVector> Phi_p);
+    ReactionPotential(std::unique_ptr<SCRF> scrf_p, std::shared_ptr<mrchem::OrbitalVector> Phi_p);
     ~ReactionPotential() override { free(NUMBER::Total); }
 
     SCRF *getHelper() { return this->helper.get(); }
diff --git a/src/qmoperators/two_electron/XCPotential.cpp b/src/qmoperators/two_electron/XCPotential.cpp
index e9be1e998..d57255f3d 100644
--- a/src/qmoperators/two_electron/XCPotential.cpp
+++ b/src/qmoperators/two_electron/XCPotential.cpp
@@ -23,9 +23,10 @@
  * <https://mrchem.readthedocs.io/>
  */
 
+#include <MRCPP/Printer>
+#include <MRCPP/Timer>
+
 #include "XCPotential.h"
-#include "MRCPP/Printer"
-#include "MRCPP/Timer"
 #include "parallel.h"
 #include "qmfunctions/Density.h"
 #include "qmfunctions/Orbital.h"
@@ -56,6 +57,11 @@ namespace mrchem {
 void XCPotential::setup(double prec) {
     if (isSetup(prec)) return;
     setApplyPrec(prec);
+    Timer timer;
+    auto plevel = Printer::getPrintLevel();
+    mrcpp::print::header(3, "Building XC operator");
+    mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
+    mrcpp::print::separator(3, '-');
     if (this->mrdft == nullptr) MSG_ERROR("XCFunctional not initialized");
     if (this->potentials.size() != 0) MSG_ERROR("Potential not properly cleared");
 
@@ -84,6 +90,19 @@ void XCPotential::setup(double prec) {
     }
 
     mrcpp::clear(xc_out, true);
+
+    if (plevel == 2) {
+        int totNodes = 0;
+        int totSize = 0;
+        for (auto i = 0; i < this->potentials.size(); i++) {
+            auto &f_i = mrcpp::get_func(this->potentials, i);
+            totNodes += f_i.getNNodes();
+            totSize += f_i.getSizeNodes();
+        }
+        auto t = timer.elapsed();
+        mrcpp::print::tree(2, "XC operator", totNodes, totSize, t);
+    }
+    mrcpp::print::footer(3, timer, 2);
 }
 
 /** @brief Clears all data in the XCPotential object */
diff --git a/src/qmoperators/two_electron/XCPotentialD1.cpp b/src/qmoperators/two_electron/XCPotentialD1.cpp
index 14cf3f895..54a58995b 100644
--- a/src/qmoperators/two_electron/XCPotentialD1.cpp
+++ b/src/qmoperators/two_electron/XCPotentialD1.cpp
@@ -58,7 +58,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Total);
             }
-            print_utils::qmfunction(2, "XC rho", rho, timer);
+            print_utils::qmfunction(3, "Compute rho", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
     } else {
@@ -70,7 +70,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Alpha);
             }
-            print_utils::qmfunction(2, "XC rho (alpha)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho (alpha)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
         { // Unperturbed beta density
@@ -81,7 +81,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Beta);
             }
-            print_utils::qmfunction(2, "XC rho (beta)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho (beta)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
     }
diff --git a/src/qmoperators/two_electron/XCPotentialD2.cpp b/src/qmoperators/two_electron/XCPotentialD2.cpp
index 62e0779f5..877f598a3 100644
--- a/src/qmoperators/two_electron/XCPotentialD2.cpp
+++ b/src/qmoperators/two_electron/XCPotentialD2.cpp
@@ -80,7 +80,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Total);
             }
-            print_utils::qmfunction(2, "XC rho_0", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_0", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
         { // Perturbed total density
@@ -91,7 +91,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Total);
             }
-            print_utils::qmfunction(2, "XC rho_1", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_1", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
     } else {
@@ -103,7 +103,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Alpha);
             }
-            print_utils::qmfunction(2, "XC rho_0 (alpha)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_0 (alpha)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
         { // Unperturbed beta density
@@ -114,7 +114,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, DensityType::Beta);
             }
-            print_utils::qmfunction(2, "XC rho_0 (beta)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_0 (beta)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
         { // Perturbed alpha density
@@ -125,7 +125,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Alpha);
             }
-            print_utils::qmfunction(2, "XC rho_1 (alpha)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_1 (alpha)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
         { // Perturbed beta density
@@ -136,7 +136,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F
                 mrcpp::copy_grid(rho.real(), grid);
                 density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Beta);
             }
-            print_utils::qmfunction(2, "XC rho_1 (beta)", rho, timer);
+            print_utils::qmfunction(3, "Compute rho_1 (beta)", rho, timer);
             dens_vec.push_back(std::make_tuple(1.0, &rho.real()));
         }
     }
diff --git a/src/scf_solver/Accelerator.cpp b/src/scf_solver/Accelerator.cpp
index 12aa165d6..25f20fdf2 100644
--- a/src/scf_solver/Accelerator.cpp
+++ b/src/scf_solver/Accelerator.cpp
@@ -107,7 +107,7 @@ void Accelerator::rotate(const ComplexMatrix &U, bool all) {
         F = U.adjoint() * F * U;
         dF = U.adjoint() * dF * U;
     }
-    mrcpp::print::time(2, "Rotating iterative subspace", t_tot);
+    mrcpp::print::time(this->pl + 2, "Rotating iterative subspace", t_tot);
 }
 
 /** @brief Update iterative history with the latest orbitals and updates
@@ -137,7 +137,7 @@ void Accelerator::push_back(OrbitalVector &Phi, OrbitalVector &dPhi, ComplexMatr
     if (historyIsFull and this->dFock.size() > 0) this->dFock.pop_front();
 
     if (not verifyOverlap(Phi)) {
-        println(2, " Clearing accelerator");
+        println(this->pl + 2, " Clearing accelerator");
         this->clear();
     }
 
@@ -146,7 +146,7 @@ void Accelerator::push_back(OrbitalVector &Phi, OrbitalVector &dPhi, ComplexMatr
     if (F != nullptr) this->fock.push_back(*F);
     if (dF != nullptr) this->dFock.push_back(*dF);
 
-    mrcpp::print::time(2, "Push back orbitals", t_tot);
+    mrcpp::print::time(this->pl + 2, "Push back orbitals", t_tot);
 }
 
 /** @brief Verify that the orbital overlap between the two last iterations is positive.
@@ -170,7 +170,7 @@ bool Accelerator::verifyOverlap(OrbitalVector &Phi) {
                 auto sqNorm = phi_i.squaredNorm();
                 auto overlap = orbital::dot(phi_i, last_i);
                 if (std::abs(overlap) < 0.5 * sqNorm) {
-                    mrcpp::print::value(2, "Overlap not verified ", std::abs(overlap));
+                    mrcpp::print::value(this->pl + 2, "Overlap not verified ", std::abs(overlap));
                     out(i) = 1;
                 }
             }
@@ -201,8 +201,8 @@ void Accelerator::accelerate(double prec, OrbitalVector &Phi, OrbitalVector &dPh
     if (this->maxHistory < 1) return;
 
     Timer t_tot;
-    auto plevel = Printer::getPrintLevel();
-    mrcpp::print::header(2, "Iterative subspace accelerator");
+    auto plevel = Printer::getPrintLevel(); // global print level
+    mrcpp::print::header(this->pl + 2, "Iterative subspace accelerator");
 
     // Deep copy into history
     this->push_back(Phi, dPhi, F, dF);
@@ -216,8 +216,8 @@ void Accelerator::accelerate(double prec, OrbitalVector &Phi, OrbitalVector &dPh
         clearLinearSystem();
     }
     printSizeNodes();
-    mrcpp::print::footer(2, t_tot, 2);
-    if (plevel == 1) mrcpp::print::time(1, "Iterative subspace accelerator", t_tot);
+    mrcpp::print::footer(this->pl + 2, t_tot, 2);
+    if (plevel == 1) mrcpp::print::time(this->pl + 1, "Iterative subspace accelerator", t_tot);
 }
 
 /** @brief Solve the linear system of equations
@@ -232,7 +232,7 @@ void Accelerator::solveLinearSystem() {
         auto tmpC = this->A[n].colPivHouseholderQr().solve(this->b[n]).eval();
         this->c.push_back(tmpC);
     }
-    mrcpp::print::time(2, "Solve linear system", t_tot);
+    mrcpp::print::time(this->pl + 2, "Solve linear system", t_tot);
 }
 
 /** @brief Copy orbitals at the given point in history into the given set.
@@ -249,7 +249,7 @@ void Accelerator::copyOrbitals(OrbitalVector &Phi, int nHistory) {
     if (nHistory >= totHistory or nHistory < 0) MSG_ABORT("Requested orbitals unavailable");
     int n = totHistory - 1 - nHistory;
     Phi = orbital::deep_copy(this->orbitals[n]);
-    mrcpp::print::time(2, "Copy orbitals", t_tot);
+    mrcpp::print::time(this->pl + 2, "Copy orbitals", t_tot);
 }
 
 /** @brief Copy orbital updates at the given point in history into the given set.
@@ -265,7 +265,7 @@ void Accelerator::copyOrbitalUpdates(OrbitalVector &dPhi, int nHistory) {
     if (nHistory >= totHistory or nHistory < 0) MSG_ABORT("Requested orbitals unavailable");
     int n = totHistory - 1 - nHistory;
     dPhi = orbital::deep_copy(this->dOrbitals[n]);
-    mrcpp::print::time(2, "Copy orbital updates", t_tot);
+    mrcpp::print::time(this->pl + 2, "Copy orbital updates", t_tot);
 }
 
 /** @brief Replaces the orbital set from a given point in history.
@@ -355,9 +355,9 @@ void Accelerator::printSizeNodes() const {
         dn += orbital::get_n_nodes(orbs_i);
         dm += orbital::get_size_nodes(orbs_i);
     }
-    mrcpp::print::separator(2, '-');
-    mrcpp::print::tree(2, "Orbital sizes", n, m, 0.0);
-    mrcpp::print::tree(2, "Update sizes", dn, dm, 0.0);
+    mrcpp::print::separator(this->pl + 2, '-');
+    mrcpp::print::tree(this->pl + 2, "Orbital sizes", n, m, 0.0);
+    mrcpp::print::tree(this->pl + 2, "Update sizes", dn, dm, 0.0);
 }
 
 } // namespace mrchem
diff --git a/src/scf_solver/Accelerator.h b/src/scf_solver/Accelerator.h
index 23ba8bafc..17d0d85eb 100644
--- a/src/scf_solver/Accelerator.h
+++ b/src/scf_solver/Accelerator.h
@@ -58,6 +58,7 @@ class Accelerator {
     virtual ~Accelerator() { this->clear(); }
     void clear();
 
+    void setLocalPrintLevel(int p) { this->pl = p; }
     void setMaxHistory(int max) { this->maxHistory = max; }
     void setMinHistory(int min) { this->minHistory = min; }
 
@@ -79,6 +80,7 @@ class Accelerator {
     void printSizeNodes() const;
 
 protected:
+    int pl{0};        ///< Local print level, relative to global Printer
     int minHistory;   ///< Accelerator is activated when history reaches this size
     int maxHistory;   ///< Oldest iteration is discarded when history exceeds this size
     bool sepOrbitals; ///< Use separate subspace for each orbital
diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp
index a211ddf05..05b7a31a8 100644
--- a/src/scf_solver/GroundStateSolver.cpp
+++ b/src/scf_solver/GroundStateSolver.cpp
@@ -171,6 +171,7 @@ void GroundStateSolver::printParameters(const std::string &calculation) const {
     print_utils::text(0, "Calculation        ", calculation);
     print_utils::text(0, "Method             ", this->methodName);
     print_utils::text(0, "Relativity         ", this->relativityName);
+    print_utils::text(0, "Environment        ", this->environmentName);
     print_utils::text(0, "Checkpointing      ", (this->checkpoint) ? "On" : "Off");
     print_utils::text(0, "Max iterations     ", o_iter.str());
     print_utils::text(0, "KAIN solver        ", o_kain.str());
@@ -217,6 +218,7 @@ void GroundStateSolver::reset() {
  */
 json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F) {
     printParameters("Optimize ground state orbitals");
+
     Timer t_tot;
     json json_out;
 
diff --git a/src/scf_solver/KAIN.cpp b/src/scf_solver/KAIN.cpp
index 24db10c8e..16d29a58b 100644
--- a/src/scf_solver/KAIN.cpp
+++ b/src/scf_solver/KAIN.cpp
@@ -120,7 +120,7 @@ void KAIN::setupLinearSystem() {
     }
 
     sortLinearSystem(A_matrices, b_vectors);
-    mrcpp::print::time(2, "Setup linear system", t_tot);
+    mrcpp::print::time(this->pl + 2, "Setup linear system", t_tot);
 }
 
 /** @brief Compute the next step for orbitals and orbital updates
@@ -213,7 +213,7 @@ void KAIN::expandSolution(double prec, OrbitalVector &Phi, OrbitalVector &dPhi,
         // F is unchanged, dF is overwritten
         *dF = fockStep;
     }
-    mrcpp::print::time(2, "Expand solution", t_tot);
+    mrcpp::print::time(this->pl + 2, "Expand solution", t_tot);
 }
 
 } // namespace mrchem
diff --git a/src/scf_solver/SCFSolver.cpp b/src/scf_solver/SCFSolver.cpp
index de2b16e4a..152129f7d 100644
--- a/src/scf_solver/SCFSolver.cpp
+++ b/src/scf_solver/SCFSolver.cpp
@@ -82,10 +82,10 @@ double SCFSolver::adjustPrecision(double error) {
     this->orbPrec[0] = std::max(this->orbPrec[0], this->orbPrec[2]);
 
     mrcpp::print::separator(2, '=');
-    mrcpp::print::value(1, "Current precision", this->orbPrec[0], "", 5);
+    mrcpp::print::value(1, "Current precision", this->orbPrec[0], "(rel)", 5);
     mrcpp::print::separator(1, '-');
-    mrcpp::print::value(2, "Orbital threshold", this->orbThrs, "", 5);
-    mrcpp::print::value(2, "Property threshold", this->propThrs, "", 5);
+    mrcpp::print::value(2, "Orbital threshold", this->orbThrs, "(abs)", 5);
+    mrcpp::print::value(2, "Property threshold", this->propThrs, "(abs)", 5);
     mrcpp::print::separator(2, '=', 2);
     return this->orbPrec[0];
 }
diff --git a/src/scf_solver/SCFSolver.h b/src/scf_solver/SCFSolver.h
index a875980fa..b6e8fe4c7 100644
--- a/src/scf_solver/SCFSolver.h
+++ b/src/scf_solver/SCFSolver.h
@@ -58,6 +58,7 @@ class SCFSolver {
     void setMaxIterations(int iter) { this->maxIter = iter; }
     void setMethodName(const std::string &name) { this->methodName = name; }
     void setRelativityName(const std::string &name) { this->relativityName = name; }
+    void setEnvironmentName(const std::string &name) { this->environmentName = name; }
 
 protected:
     int history{0};                      ///< Maximum length of KAIN history
@@ -69,6 +70,7 @@ class SCFSolver {
     double orbPrec[3]{-1.0, -1.0, -1.0}; ///< Dynamic precision: [current_prec, start_prec, end_prec]
     std::string methodName;              ///< Name of electronic structure method to appear in output
     std::string relativityName{"None"};  ///< Name of ZORA method
+    std::string environmentName{"None"}; ///< Aggregated name for external environment
 
     std::vector<double> error;    ///< Convergence orbital error
     std::vector<double> property; ///< Convergence property error
diff --git a/src/tensor/RankZeroOperator.cpp b/src/tensor/RankZeroOperator.cpp
index e1f89c099..656c3a7b3 100644
--- a/src/tensor/RankZeroOperator.cpp
+++ b/src/tensor/RankZeroOperator.cpp
@@ -292,7 +292,7 @@ OrbitalVector RankZeroOperator::operator()(OrbitalVector &inp) {
         out.push_back(out_i);
         std::stringstream o_name;
         o_name << O.name() << "|" << i << ">";
-        print_utils::qmfunction(3, o_name.str(), out_i, t1);
+        print_utils::qmfunction(4, o_name.str(), out_i, t1);
     }
     return out;
 }
@@ -312,7 +312,7 @@ OrbitalVector RankZeroOperator::dagger(OrbitalVector &inp) {
         out.push_back(out_i);
         std::stringstream o_name;
         o_name << O.name() << "^dagger|" << i << ">";
-        print_utils::qmfunction(3, o_name.str(), out_i, t1);
+        print_utils::qmfunction(4, o_name.str(), out_i, t1);
     }
     return out;
 }
diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp
index d2bb985fb..65e4e6626 100644
--- a/src/utils/print_utils.cpp
+++ b/src/utils/print_utils.cpp
@@ -72,6 +72,13 @@ std::string print_utils::dbl_to_str(double d, int p, bool sci) {
     }
     return o_dbl.str();
 }
+/** @brief Prints text centered according to current print width */
+void print_utils::centered_text(int level, const std::string &txt) {
+    int pwidth = Printer::getWidth();
+    int txt_width = txt.size();
+    int pre_spaces = (pwidth - txt_width) / 2;
+    println(level, std::string(pre_spaces, ' ') << txt);
+}
 
 void print_utils::headline(int level, const std::string &txt) {
     auto pwidth = Printer::getWidth();
@@ -91,18 +98,23 @@ void print_utils::headline(int level, const std::string &txt) {
     mrcpp::print::separator(level, ' ');
 }
 
-void print_utils::text(int level, const std::string &txt, const std::string &val) {
+void print_utils::text(int level, const std::string &txt, const std::string &val, bool ralign) {
     int w0 = Printer::getWidth() - 2;
     int w1 = w0 * 2 / 9;
     int w2 = w0 - 3 * w1;
     int w3 = w2 - (txt.size() + 1);
 
+    // Right-aligning
+    int shift;
+    shift = (ralign) ? w0 - w2 - val.size() - 1 : 0;
+
     std::stringstream o;
-    o << " " << txt << std::string(w3, ' ') << ": " << val;
+    o << " " << txt << std::string(w3, ' ') << ": " << std::string(shift, ' ') << val;
     println(level, o.str());
 }
 
 void print_utils::json(int level, const nlohmann::json &j, bool ralign) {
+    if (level > mrcpp::Printer::getPrintLevel()) return;
     // Determine longest name
     // This will be used for the spacing left of :
     // if the name is too long to be aligned with
@@ -127,7 +139,7 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) {
         int w0 = Printer::getWidth() - 2; // Defines the printable width
         int w1 = w0 * 2 / 9;              // Space dedicated to the json key
         int w2 = w0 - 3 * w1;
-        int w3 = w2 - (val.size() + 1);
+        // int w3 = w2 - (val.size() + 1);
 
         // Some paddings for book keeping
         int frontEndPadding = 2; // Empty space at beginning and end
diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h
index 25127a510..2b29cc020 100644
--- a/src/utils/print_utils.h
+++ b/src/utils/print_utils.h
@@ -33,8 +33,9 @@
 namespace mrchem {
 namespace print_utils {
 
+void centered_text(int level, const std::string &txt);
 void headline(int level, const std::string &txt);
-void text(int level, const std::string &txt, const std::string &val);
+void text(int level, const std::string &txt, const std::string &val, bool ralign = false);
 void json(int level, const nlohmann::json &j, bool ralign = false);
 void coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p = -1, bool s = false);
 void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false);
diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp
index 60e553d77..23dc9d234 100644
--- a/tests/solventeffect/reaction_operator.cpp
+++ b/tests/solventeffect/reaction_operator.cpp
@@ -90,7 +90,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") {
     int kain = 4;
     auto scrf_p = std::make_unique<SCRF>(dielectric_func, molecule, P_p, D_p, prec, kain, 100, true, "dynamic", "scrf", "total");
 
-    auto Reo = std::make_shared<ReactionOperator>(scrf_p, Phi_p);
+    auto Reo = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
     Reo->setTesting();
     Reo->setup(prec);
     double total_energy = Reo->getTotalEnergy();