diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst
index d4775f9b2..ed857ac12 100644
--- a/doc/users/user_ref.rst
+++ b/doc/users/user_ref.rst
@@ -845,15 +845,6 @@ User input reference
       
         **Default** ``True``
       
-       :optimizer: Choose which function to use in the KAIN solver, the surface charge ``density`` (gamma) or the reaction ``potential`` (V_R). 
-      
-        **Type** ``str``
-      
-        **Default** ``potential``
-      
-        **Predicates**
-          - ``value.lower() in ['density', 'potential']``
-      
        :density_type: What part of the total molecular charge density to use in the algorithm. ``total`` uses the total charge density. ``nuclear`` uses only the nuclear part of the total charge density. ``electronic`` uses only the electronic part of the total charge density. 
       
         **Type** ``str``
diff --git a/external/upstream/fetch_mrcpp.cmake b/external/upstream/fetch_mrcpp.cmake
index c70eb4db9..d74c6d193 100644
--- a/external/upstream/fetch_mrcpp.cmake
+++ b/external/upstream/fetch_mrcpp.cmake
@@ -39,7 +39,7 @@ else()
     GIT_REPOSITORY
       https://github.com/MRChemSoft/mrcpp.git
     GIT_TAG
-      a3618d1498410124ec47bef777adf8b50a7e71b0
+      f8def0a086da6410e5dd8e078de4f6b6305b6ea3
   )
 
   FetchContent_GetProperties(mrcpp_sources)
diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py
index 178d38271..506c4f52b 100644
--- a/python/mrchem/helpers.py
+++ b/python/mrchem/helpers.py
@@ -77,7 +77,6 @@ def write_scf_fock(user_dict, wf_dict, origin):
             "poisson_prec": user_dict["world_prec"],
             "kain": user_dict["PCM"]["SCRF"]["kain"],
             "max_iter": user_dict["PCM"]["SCRF"]["max_iter"],
-            "optimizer": user_dict["PCM"]["SCRF"]["optimizer"],
             "dynamic_thrs": user_dict["PCM"]["SCRF"]["dynamic_thrs"],
             "density_type": user_dict["PCM"]["SCRF"]["density_type"],
             "epsilon_in": user_dict["PCM"]["Permittivity"]["epsilon_in"],
diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py
index 22657aab0..fe1fb868b 100644
--- a/python/mrchem/input_parser/api.py
+++ b/python/mrchem/input_parser/api.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# This file was automatically generated by parselglossy on 2023-09-25
+# This file was automatically generated by parselglossy on 2023-10-30
 # Editing is *STRONGLY DISCOURAGED*
 
 from copy import deepcopy
@@ -572,13 +572,6 @@ def stencil() -> JSONDict:
                                                             {   'default': True,
                                                                 'name': 'dynamic_thrs',
                                                                 'type': 'bool'},
-                                                            {   'default': 'potential',
-                                                                'name': 'optimizer',
-                                                                'predicates': [   'value.lower() '
-                                                                                  'in '
-                                                                                  "['density', "
-                                                                                  "'potential']"],
-                                                                'type': 'str'},
                                                             {   'default': 'total',
                                                                 'name': 'density_type',
                                                                 'predicates': [   'value.lower() '
diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst
index d4775f9b2..ed857ac12 100644
--- a/python/mrchem/input_parser/docs/user_ref.rst
+++ b/python/mrchem/input_parser/docs/user_ref.rst
@@ -845,15 +845,6 @@ User input reference
       
         **Default** ``True``
       
-       :optimizer: Choose which function to use in the KAIN solver, the surface charge ``density`` (gamma) or the reaction ``potential`` (V_R). 
-      
-        **Type** ``str``
-      
-        **Default** ``potential``
-      
-        **Predicates**
-          - ``value.lower() in ['density', 'potential']``
-      
        :density_type: What part of the total molecular charge density to use in the algorithm. ``total`` uses the total charge density. ``nuclear`` uses only the nuclear part of the total charge density. ``electronic`` uses only the electronic part of the total charge density. 
       
         **Type** ``str``
diff --git a/python/template.yml b/python/template.yml
index 36740f7ee..09b3021dc 100644
--- a/python/template.yml
+++ b/python/template.yml
@@ -828,14 +828,6 @@ sections:
               are close to convergence (``mo_residual < world_prec*10``) the convergence
               threshold will be set equal to ``world_prec``.
               ``false`` uses ``world_prec`` as convergence threshold throughout.
-          - name: optimizer
-            type: str
-            default: potential
-            predicates:
-              - value.lower() in ['density', 'potential']
-            docstring: |
-              Choose which function to use in the KAIN solver, the surface charge ``density``
-              (gamma) or the reaction ``potential`` (V_R).
           - name: kain
             type: int
             default: user['SCF']['kain']
diff --git a/src/driver.cpp b/src/driver.cpp
index 149062fac..cdccf7fc0 100644
--- a/src/driver.cpp
+++ b/src/driver.cpp
@@ -44,6 +44,7 @@
 #include "utils/math_utils.h"
 #include "utils/print_utils.h"
 
+#include "chemistry/chemistry_utils.h"
 #include "qmfunctions/Density.h"
 #include "qmfunctions/Orbital.h"
 #include "qmfunctions/density_utils.h"
@@ -1050,18 +1051,18 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
 
         auto kain = json_fock["reaction_operator"]["kain"];
         auto max_iter = json_fock["reaction_operator"]["max_iter"];
-        auto optimizer = json_fock["reaction_operator"]["optimizer"];
         auto dynamic_thrs = json_fock["reaction_operator"]["dynamic_thrs"];
         auto density_type = json_fock["reaction_operator"]["density_type"];
         auto eps_i = json_fock["reaction_operator"]["epsilon_in"];
         auto eps_o = json_fock["reaction_operator"]["epsilon_out"];
         auto formulation = json_fock["reaction_operator"]["formulation"];
-        auto accelerate_pot = (optimizer == "potential") ? true : false;
 
         Permittivity dielectric_func(*cavity_p, eps_i, eps_o, formulation);
         dielectric_func.printParameters();
+        Density rho_nuc(false);
+        rho_nuc = chemistry::compute_nuclear_density(poisson_prec, nuclei, 100);
 
-        auto scrf_p = std::make_unique<SCRF>(dielectric_func, nuclei, P_p, D_p, poisson_prec, kain, max_iter, accelerate_pot, dynamic_thrs, density_type);
+        auto scrf_p = std::make_unique<SCRF>(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type);
         auto V_R = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
         F.getReactionOperator() = V_R;
     }
diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp
index 3fbfcb902..7ac30fc6e 100644
--- a/src/environment/SCRF.cpp
+++ b/src/environment/SCRF.cpp
@@ -30,8 +30,6 @@
 #include <MRCPP/Printer>
 #include <MRCPP/Timer>
 
-#include "chemistry/PhysicalConstants.h"
-#include "chemistry/chemistry_utils.h"
 #include "qmfunctions/density_utils.h"
 #include "qmoperators/two_electron/ReactionPotential.h"
 #include "scf_solver/KAIN.h"
@@ -48,37 +46,24 @@ using OrbitalVector_p = std::shared_ptr<mrchem::OrbitalVector>;
 
 namespace mrchem {
 
-SCRF::SCRF(Permittivity e, const Nuclei &N, PoissonOperator_p P, DerivativeOperator_p D, double orb_prec, int kain_hist, int max_iter, bool acc_pot, bool dyn_thrs, std::string density_type)
-        : accelerate_Vr(acc_pot)
-        , dynamic_thrs(dyn_thrs)
+SCRF::SCRF(const Permittivity &e, const Density &rho_nuc, PoissonOperator_p P, DerivativeOperator_p D, int kain_hist, int max_iter, bool dyn_thrs, const std::string &density_type)
+        : dynamic_thrs(dyn_thrs)
         , density_type(density_type)
         , max_iter(max_iter)
         , history(kain_hist)
-        , apply_prec(orb_prec)
-        , conv_thrs(1.0)
-        , mo_residual(1.0)
         , epsilon(e)
-        , rho_nuc(false)
-        , rho_ext(false)
-        , rho_tot(false)
+        , rho_nuc(rho_nuc)
         , Vr_n(false)
-        , dVr_n(false)
-        , Vr_nm1(false)
-        , gamma_n(false)
-        , dgamma_n(false)
-        , gamma_nm1(false)
         , derivative(D)
-        , poisson(P) {
-    setDCavity();
-    rho_nuc = chemistry::compute_nuclear_density(this->apply_prec, N, 1000);
-}
+        , poisson(P) {}
 
 SCRF::~SCRF() {
-    mrcpp::clear(this->d_cavity, true);
+    this->rho_nuc.free(NUMBER::Real);
+    clear();
 }
 
 void SCRF::clear() {
-    this->rho_tot.free(NUMBER::Real);
+    this->apply_prec = -1.0;
 }
 
 double SCRF::setConvergenceThreshold(double prec) {
@@ -88,59 +73,65 @@ double SCRF::setConvergenceThreshold(double prec) {
     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);
-    mrcpp::FunctionTree<3> *dz_cavity = new mrcpp::FunctionTree<3>(*MRA);
-    d_cavity.push_back(std::make_tuple(1.0, dx_cavity));
-    d_cavity.push_back(std::make_tuple(1.0, dy_cavity));
-    d_cavity.push_back(std::make_tuple(1.0, dz_cavity));
-    mrcpp::project<3>(this->apply_prec / 100, this->d_cavity, this->epsilon.getGradVector());
-}
-
-void SCRF::computeDensities(OrbitalVector &Phi) {
+void SCRF::computeDensities(OrbitalVector &Phi, Density &rho_out) {
     Timer timer;
-    resetComplexFunction(this->rho_tot);
+
     Density rho_el(false);
     density::compute(this->apply_prec, rho_el, Phi, DensityType::Total);
     rho_el.rescale(-1.0);
+
     if (this->density_type == "electronic") {
-        mrcpp::cplxfunc::deep_copy(this->rho_tot, rho_el);
+        mrcpp::cplxfunc::deep_copy(rho_out, rho_el);
     } else if (this->density_type == "nuclear") {
-        mrcpp::cplxfunc::deep_copy(this->rho_tot, this->rho_nuc);
+        mrcpp::cplxfunc::deep_copy(rho_out, this->rho_nuc);
     } else {
-        mrcpp::cplxfunc::add(this->rho_tot, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); // probably change this into a vector
+        mrcpp::cplxfunc::add(rho_out, 1.0, rho_el, 1.0, this->rho_nuc, -1.0);
     }
-    print_utils::qmfunction(3, "Vacuum density", this->rho_tot, timer);
+    print_utils::qmfunction(3, "Vacuum density", rho_out, timer);
 }
 
 void SCRF::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) {
-    auto d_V = mrcpp::gradient(*derivative, potential.real());
+
+    auto d_V = mrcpp::gradient(*derivative, potential.real()); // FunctionTreeVector
     resetComplexFunction(out_gamma);
-    mrcpp::dot(this->apply_prec, out_gamma.real(), d_V, this->d_cavity);
+
+    for (int d = 0; d < 3; d++) {
+        mrcpp::AnalyticFunction<3> d_cav(this->epsilon.getGradVector()[d]);
+        mrcpp::ComplexFunction cplxfunc_prod;
+        mrcpp::cplxfunc::multiply(cplxfunc_prod, get_func(d_V, d), d_cav, this->apply_prec, 1);
+        // add result into out_gamma
+        if (d == 0) {
+            mrcpp::cplxfunc::deep_copy(out_gamma, cplxfunc_prod);
+        } else {
+            out_gamma.add(1.0, cplxfunc_prod);
+        }
+    }
+
     out_gamma.rescale(std::log((epsilon.getEpsIn() / epsilon.getEpsOut())) * (1.0 / (4.0 * mrcpp::pi)));
     mrcpp::clear(d_V, true);
 }
 
-mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma) {
+mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, mrchem::OrbitalVector &Phi) {
     mrcpp::ComplexFunction Poisson_func;
     mrcpp::ComplexFunction rho_eff;
     mrcpp::ComplexFunction first_term;
-    mrcpp::ComplexFunction Vr;
-    mrcpp::ComplexFunction eps_inv;
-    eps_inv.alloc(NUMBER::Real);
-    Vr.alloc(NUMBER::Real);
+    mrcpp::ComplexFunction Vr_np1;
+    Vr_np1.alloc(NUMBER::Real);
 
     this->epsilon.flipFunction(true);
-    mrcpp::cplxfunc::project(eps_inv, this->epsilon, NUMBER::Real, this->apply_prec / 100);
+
+    Density rho_tot(false);
+    computeDensities(Phi, rho_tot);
+
+    mrcpp::cplxfunc::multiply(first_term, rho_tot, this->epsilon, this->apply_prec);
     this->epsilon.flipFunction(false);
-    mrcpp::cplxfunc::multiply(first_term, this->rho_tot, eps_inv, this->apply_prec);
-    mrcpp::cplxfunc::add(rho_eff, 1.0, first_term, -1.0, this->rho_tot, -1.0);
 
-    mrcpp::cplxfunc::add(Poisson_func, 1.0, in_gamma, 1.0, rho_eff, -1.0);
-    mrcpp::apply(this->apply_prec, Vr.real(), *poisson, Poisson_func.real());
+    mrcpp::cplxfunc::add(rho_eff, 1.0, first_term, -1.0, rho_tot, -1.0);
+    rho_tot.free(NUMBER::Real);
 
-    return Vr;
+    mrcpp::cplxfunc::add(Poisson_func, 1.0, in_gamma, 1.0, rho_eff, -1.0);
+    mrcpp::apply(this->apply_prec, Vr_np1.real(), *poisson, Poisson_func.real());
+    return Vr_np1;
 }
 
 void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain) {
@@ -152,6 +143,9 @@ void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFu
     phi_n[0] = func;
     dPhi_n[0] = dfunc;
 
+    phi_n[0].setRank(mrcpp::mpi::wrk_rank);
+    dPhi_n[0].setRank(mrcpp::mpi::wrk_rank);
+
     kain.accelerate(this->apply_prec, phi_n, dPhi_n);
 
     func = phi_n[0];
@@ -161,7 +155,7 @@ void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFu
     dPhi_n.clear();
 }
 
-void SCRF::nestedSCRF(mrcpp::ComplexFunction V_vac) {
+void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, std::shared_ptr<mrchem::OrbitalVector> Phi_p) {
     KAIN kain(this->history);
     kain.setLocalPrintLevel(10);
 
@@ -172,47 +166,37 @@ void SCRF::nestedSCRF(mrcpp::ComplexFunction V_vac) {
     while (update >= this->conv_thrs && iter <= max_iter) {
         Timer t_iter;
         // solve the poisson equation
-        mrcpp::ComplexFunction Vr_np1 = solvePoissonEquation(this->gamma_n);
+        mrcpp::ComplexFunction V_tot;
+        mrcpp::ComplexFunction gamma_n;
+        mrcpp::ComplexFunction dVr_n;
+
+        mrcpp::cplxfunc::add(V_tot, 1.0, this->Vr_n, 1.0, V_vac, -1.0);
+        computeGamma(V_tot, gamma_n);
+        auto Vr_np1 = solvePoissonEquation(gamma_n, *Phi_p);
         norm = Vr_np1.norm();
 
         // use a convergence accelerator
-        resetComplexFunction(this->dVr_n);
-        mrcpp::cplxfunc::add(this->dVr_n, 1.0, Vr_np1, -1.0, this->Vr_n, -1.0);
+
+        mrcpp::cplxfunc::add(dVr_n, 1.0, Vr_np1, -1.0, this->Vr_n, -1.0);
         update = dVr_n.norm();
 
-        if (iter > 1 and this->history > 0 and this->accelerate_Vr) {
+        if (iter > 1 and this->history > 0) {
             accelerateConvergence(dVr_n, Vr_n, kain);
             Vr_np1.free(NUMBER::Real);
             mrcpp::cplxfunc::add(Vr_np1, 1.0, Vr_n, 1.0, dVr_n, -1.0);
         }
 
         // set up for next iteration
-        mrcpp::ComplexFunction V_tot;
-        mrcpp::cplxfunc::add(V_tot, 1.0, Vr_np1, 1.0, V_vac, -1.0);
-        updateCurrentReactionPotential(Vr_np1); // push_back() maybe
-
-        mrcpp::ComplexFunction gamma_np1;
-        computeGamma(V_tot, gamma_np1);
-
-        resetComplexFunction(dgamma_n);
-        mrcpp::cplxfunc::add(this->dgamma_n, 1.0, gamma_np1, -1.0, this->gamma_n, -1.0);
-
-        if (iter > 1 and this->history > 0 and (not this->accelerate_Vr)) {
-            accelerateConvergence(dgamma_n, gamma_n, kain);
-            gamma_np1.free(NUMBER::Real);
-            mrcpp::cplxfunc::add(gamma_np1, 1.0, gamma_n, 1.0, dgamma_n, -1.0);
-        }
+        resetComplexFunction(this->Vr_n);
+        mrcpp::cplxfunc::deep_copy(this->Vr_n, Vr_np1);
+        Vr_np1.free(NUMBER::Real);
 
-        updateCurrentGamma(gamma_np1);
         printConvergenceRow(iter, norm, update, t_iter.elapsed());
         iter++;
     }
     if (iter > max_iter) println(0, "Reaction potential failed to converge after " << iter - 1 << " iterations, residual " << update);
     mrcpp::print::separator(3, '-');
-    this->dVr_n.real().clear();
-    this->dVr_n.real().setZero();
-    this->dgamma_n.real().clear();
-    this->dgamma_n.real().setZero();
+
     kain.clear();
 }
 
@@ -242,63 +226,38 @@ void SCRF::printConvergenceRow(int i, double norm, double update, double time) c
 
 mrcpp::ComplexFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) {
     this->apply_prec = prec;
-    computeDensities(*Phi);
+    Density rho_tot(false);
+    computeDensities(*Phi, rho_tot);
     Timer t_vac;
     mrcpp::ComplexFunction V_vac;
     V_vac.alloc(NUMBER::Real);
-    mrcpp::apply(this->apply_prec, V_vac.real(), *poisson, this->rho_tot.real());
+    mrcpp::apply(this->apply_prec, V_vac.real(), *poisson, rho_tot.real());
+    rho_tot.free(NUMBER::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())) {
+    if (not this->Vr_n.hasReal()) {
         mrcpp::ComplexFunction gamma_0;
         mrcpp::ComplexFunction V_tot;
         computeGamma(V_vac, gamma_0);
-        this->Vr_n = solvePoissonEquation(gamma_0);
-        mrcpp::cplxfunc::add(V_tot, 1.0, V_vac, 1.0, this->Vr_n, -1.0);
-        computeGamma(V_tot, this->gamma_n);
+        this->Vr_n = solvePoissonEquation(gamma_0, *Phi);
     }
 
     // update the potential/gamma before doing anything with them
 
-    if (accelerate_Vr) {
-        mrcpp::ComplexFunction temp_Vr_n;
-        mrcpp::cplxfunc::add(temp_Vr_n, 1.0, this->Vr_n, 1.0, this->dVr_n, -1.0);
-        mrcpp::cplxfunc::deep_copy(this->Vr_n, temp_Vr_n);
-        temp_Vr_n.free(NUMBER::Real);
-        mrcpp::ComplexFunction V_tot;
-        mrcpp::cplxfunc::add(V_tot, 1.0, this->Vr_n, 1.0, V_vac, -1.0);
-        resetComplexFunction(this->gamma_n);
-        computeGamma(V_tot, this->gamma_n);
-    } else {
-        mrcpp::ComplexFunction temp_gamma_n;
-        mrcpp::cplxfunc::add(temp_gamma_n, 1.0, this->gamma_n, 1.0, this->dgamma_n, -1.0);
-        mrcpp::cplxfunc::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;
-    nestedSCRF(V_vac);
+    nestedSCRF(V_vac, Phi);
     print_utils::qmfunction(3, "Reaction potential", this->Vr_n, t_scrf);
     return this->Vr_n;
 }
 
-double SCRF::getNuclearEnergy() {
-    return mrcpp::cplxfunc::dot(this->rho_nuc, this->Vr_n).real();
-}
-
-double SCRF::getElectronicEnergy() {
-    mrcpp::ComplexFunction rho_el;
-    rho_el.alloc(NUMBER::Real);
-    mrcpp::cplxfunc::add(rho_el, 1.0, this->rho_tot, -1.0, this->rho_nuc, -1.0);
-    return mrcpp::cplxfunc::dot(rho_el, this->Vr_n).real();
-}
+auto SCRF::computeEnergies(const Density &rho_el) -> std::tuple<double, double> {
+    auto Er_nuc = 0.5 * mrcpp::cplxfunc::dot(this->rho_nuc, this->Vr_n).real();
 
-double SCRF::getTotalEnergy() {
-    return mrcpp::cplxfunc::dot(this->rho_tot, this->Vr_n).real();
+    auto Er_el = 0.5 * mrcpp::cplxfunc::dot(rho_el, this->Vr_n).real();
+    return std::make_tuple(Er_el, Er_nuc);
 }
 
 void SCRF::resetComplexFunction(mrcpp::ComplexFunction &function) {
@@ -307,22 +266,6 @@ void SCRF::resetComplexFunction(mrcpp::ComplexFunction &function) {
     function.alloc(NUMBER::Real);
 }
 
-void SCRF::updateCurrentReactionPotential(mrcpp::ComplexFunction &Vr_np1) {
-    resetComplexFunction(this->Vr_nm1);
-    mrcpp::cplxfunc::deep_copy(this->Vr_nm1, this->Vr_n);
-    resetComplexFunction(this->Vr_n);
-    mrcpp::cplxfunc::deep_copy(this->Vr_n, Vr_np1);
-    Vr_np1.free(NUMBER::Real);
-}
-
-void SCRF::updateCurrentGamma(mrcpp::ComplexFunction &gamma_np1) {
-    resetComplexFunction(this->gamma_nm1);
-    mrcpp::cplxfunc::deep_copy(this->gamma_nm1, this->gamma_n);
-    resetComplexFunction(this->gamma_n);
-    mrcpp::cplxfunc::deep_copy(this->gamma_n, gamma_np1);
-    gamma_np1.free(NUMBER::Real);
-}
-
 void SCRF::printParameters() const {
     std::stringstream o_iter;
     if (this->max_iter > 0) {
@@ -340,7 +283,7 @@ void SCRF::printParameters() const {
 
     nlohmann::json data = {
         {"Method                ", "SCRF"},
-        {"Optimizer             ", (accelerate_Vr) ? "Potential" : "Density"},
+        {"Optimizer             ", "Potential"},
         {"Max iterations        ", max_iter},
         {"KAIN solver           ", o_kain.str()},
         {"Density type          ", density_type},
diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h
index 0a802113e..df9312841 100644
--- a/src/environment/SCRF.h
+++ b/src/environment/SCRF.h
@@ -30,93 +30,68 @@
 #include "qmfunctions/Orbital.h"
 
 namespace mrchem {
+class KAIN;
 /** @class SCRF
  *  @brief class that performs the computation of the  ReactionPotential, named Self Consistent Reaction Field.
  */
-class Nuclei;
-class KAIN;
 class SCRF final {
 public:
-    SCRF(Permittivity e,
-         const Nuclei &N,
+    SCRF(const Permittivity &e,
+         const Density &rho_nuc,
          std::shared_ptr<mrcpp::PoissonOperator> P,
          std::shared_ptr<mrcpp::DerivativeOperator<3>> D,
-         double orb_prec,
          int kain_hist,
          int max_iter,
-         bool acc_pot,
          bool dyn_thrs,
-         std::string density_type);
+         const std::string &density_type);
     ~SCRF();
-    void UpdateExternalDensity(Density new_density) { this->rho_ext = new_density; }
 
     double setConvergenceThreshold(double prec);
 
     mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->Vr_n; }
-    mrcpp::ComplexFunction &getPreviousReactionPotential() { return this->Vr_nm1; }
-    mrcpp::ComplexFunction &getCurrentDifferenceReactionPotential() { return this->dVr_n; }
-
-    mrcpp::ComplexFunction &getCurrentGamma() { return this->gamma_n; }
-    mrcpp::ComplexFunction &getPreviousGamma() { return this->gamma_nm1; }
-    mrcpp::ComplexFunction &getCurrentDifferenceGamma() { return this->dgamma_n; }
 
     Permittivity &getPermittivity() { return this->epsilon; }
 
     void updateMOResidual(double const err_t) { this->mo_residual = err_t; }
 
     friend class ReactionPotential;
+    auto computeEnergies(const Density &rho_el) -> std::tuple<double, double>;
 
 protected:
     void clear();
 
 private:
-    bool accelerate_Vr;
     bool dynamic_thrs;
     std::string density_type;
 
     int max_iter;
     int history;
-    double apply_prec;
-    double conv_thrs;
-    double mo_residual;
+    double apply_prec{-1.0};
+    double conv_thrs{1.0};
+    double mo_residual{1.0};
 
     Permittivity epsilon;
 
-    Density rho_nuc;
-    Density rho_ext;
-    Density rho_tot;
+    Density rho_nuc; // As of right now, this is the biggest memory hog.
+    // Alternative could be to precompute its contributions, as a potential is not as heavy as a density (maybe)
+    // another one could be to define a representable function which only has the exact analytical form of the nuclear contribution.
 
     mrcpp::ComplexFunction Vr_n;
-    mrcpp::ComplexFunction dVr_n;
-    mrcpp::ComplexFunction Vr_nm1;
 
-    mrcpp::ComplexFunction gamma_n;
-    mrcpp::ComplexFunction dgamma_n;
-    mrcpp::ComplexFunction gamma_nm1;
-
-    mrcpp::FunctionTreeVector<3> d_cavity; //!< Vector containing the 3 partial derivatives of the cavity function
     std::shared_ptr<mrcpp::DerivativeOperator<3>> derivative;
     std::shared_ptr<mrcpp::PoissonOperator> poisson;
 
-    void setDCavity();
-
-    void computeDensities(OrbitalVector &Phi);
+    void computeDensities(OrbitalVector &Phi, Density &rho_out);
     void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma);
 
-    mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma);
+    mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, mrchem::OrbitalVector &Phi);
 
     void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain);
 
-    // TODO    void variationalSCRF(mrcpp::ComplexFunction V_vac);
-    void nestedSCRF(mrcpp::ComplexFunction V_vac);
-    mrcpp::ComplexFunction &setup(double prec, const std::shared_ptr<mrchem::OrbitalVector> &Phi);
+    void nestedSCRF(const mrcpp::ComplexFunction &V_vac, std::shared_ptr<mrchem::OrbitalVector> Phi_p);
+    mrcpp::ComplexFunction &setup(double prec, const std::shared_ptr<mrchem::OrbitalVector> &Phi_p);
 
-    double getNuclearEnergy();
-    double getElectronicEnergy();
-    double getTotalEnergy();
     void resetComplexFunction(mrcpp::ComplexFunction &function);
-    void updateCurrentReactionPotential(mrcpp::ComplexFunction &Vr_np1);
-    void updateCurrentGamma(mrcpp::ComplexFunction &gamma_np1);
 
     void printParameters() const;
     void printConvergenceRow(int i, double norm, double update, double time) const;
diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp
index 14defd2a1..8a5a512e1 100644
--- a/src/qmoperators/two_electron/FockBuilder.cpp
+++ b/src/qmoperators/two_electron/FockBuilder.cpp
@@ -36,6 +36,7 @@
 #include "chemistry/chemistry_utils.h"
 #include "properties/SCFEnergy.h"
 #include "qmfunctions/Orbital.h"
+#include "qmfunctions/density_utils.h"
 #include "qmfunctions/orbital_utils.h"
 #include "qmoperators/one_electron/ElectricFieldOperator.h"
 #include "qmoperators/one_electron/IdentityOperator.h"
@@ -169,9 +170,15 @@ SCFEnergy FockBuilder::trace(OrbitalVector &Phi, const Nuclei &nucs) {
 
     // Reaction potential part
     if (this->Ro != nullptr) {
-        Er_nuc = 0.5 * this->Ro->getNuclearEnergy();
-        Er_tot = 0.5 * this->Ro->getTotalEnergy();
-        Er_el = 0.5 * this->Ro->getElectronicEnergy();
+        Density rho_el(false);
+        density::compute(this->prec, rho_el, Phi, DensityType::Total);
+        rho_el.rescale(-1.0);
+        std::tuple<double, double> Er = this->Ro->getHelper()->computeEnergies(rho_el);
+
+        Er_el = std::get<0>(Er);
+        Er_nuc = std::get<1>(Er);
+
+        Er_tot = Er_nuc + Er_el;
     }
 
     // Kinetic part
@@ -182,9 +189,7 @@ SCFEnergy FockBuilder::trace(OrbitalVector &Phi, const Nuclei &nucs) {
     }
 
     // Electronic part
-    if (this->nuc != nullptr) {
-        E_en = this->nuc->trace(Phi).real();
-    }
+    if (this->nuc != nullptr) { E_en = this->nuc->trace(Phi).real(); }
 
     if (this->coul != nullptr) E_ee = 0.5 * this->coul->trace(Phi).real();
     if (this->ex != nullptr) E_x = -this->exact_exchange * this->ex->trace(Phi).real();
diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h
index 53c1106f2..d5ab309cc 100644
--- a/src/qmoperators/two_electron/ReactionOperator.h
+++ b/src/qmoperators/two_electron/ReactionOperator.h
@@ -52,21 +52,11 @@ class ReactionOperator final : public RankZeroOperator {
 
     ComplexDouble trace(OrbitalVector &Phi) { return RankZeroOperator::trace(Phi); }
 
-    double getTotalEnergy() { return this->potential->getTotalEnergy(); }
-    double getNuclearEnergy() { return this->potential->getNuclearEnergy(); }
-    double getElectronicEnergy() { return this->potential->getElectronicEnergy(); }
-
     SCRF *getHelper() { return this->potential->getHelper(); }
     std::shared_ptr<ReactionPotential> getPotential() { return this->potential; }
     void updateMOResidual(double const err_t) { this->potential->updateMOResidual(err_t); }
 
     mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->potential->getCurrentReactionPotential(); }
-    mrcpp::ComplexFunction &getPreviousReactionPotential() { return this->potential->getPreviousReactionPotential(); }
-    mrcpp::ComplexFunction &getCurrentDifferenceReactionPotential() { return this->potential->getCurrentDifferenceReactionPotential(); }
-
-    mrcpp::ComplexFunction &getCurrentGamma() { return this->potential->getCurrentGamma(); }
-    mrcpp::ComplexFunction &getPreviousGamma() { return this->potential->getPreviousGamma(); }
-    mrcpp::ComplexFunction &getCurrentDifferenceGamma() { return this->potential->getCurrentDifferenceGamma(); }
 
 private:
     std::shared_ptr<ReactionPotential> potential{nullptr};
diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h
index e528401be..698fc123c 100644
--- a/src/qmoperators/two_electron/ReactionPotential.h
+++ b/src/qmoperators/two_electron/ReactionPotential.h
@@ -48,21 +48,12 @@ class ReactionPotential final : public QMPotential {
     ~ReactionPotential() override { free(NUMBER::Total); }
 
     SCRF *getHelper() { return this->helper.get(); }
-    double getElectronicEnergy() { return this->helper->getElectronicEnergy(); }
-    double getNuclearEnergy() { return this->helper->getNuclearEnergy(); }
-    double getTotalEnergy() { return this->helper->getTotalEnergy(); }
 
     /** @brief Updates the helper.mo_residual member variable. This variable is used to set the convergence criterion in
      * the dynamic convergence method. */
     void updateMOResidual(double const err_t) { this->helper->mo_residual = err_t; }
 
     mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->helper->getCurrentReactionPotential(); }
-    mrcpp::ComplexFunction &getPreviousReactionPotential() { return this->helper->getPreviousReactionPotential(); }
-    mrcpp::ComplexFunction &getCurrentDifferenceReactionPotential() { return this->helper->getCurrentDifferenceReactionPotential(); }
-
-    mrcpp::ComplexFunction &getCurrentGamma() { return this->helper->getCurrentGamma(); }
-    mrcpp::ComplexFunction &getPreviousGamma() { return this->helper->getPreviousGamma(); }
-    mrcpp::ComplexFunction &getCurrentDifferenceGamma() { return this->helper->getCurrentDifferenceGamma(); }
 
     friend class ReactionOperator;
 
diff --git a/tests/li_solv/li.inp b/tests/li_solv/li.inp
index 783f7d855..fcc9f4217 100644
--- a/tests/li_solv/li.inp
+++ b/tests/li_solv/li.inp
@@ -16,8 +16,7 @@
   "SCRF": {
     "kain": 5,
     "max_iter": 100,
-    "dynamic_thrs": false,
-    "optimizer": "potential"
+    "dynamic_thrs": false
   },
   "Cavity": {
     "spheres": "0 4.0 1.0 0.0 0.5"
diff --git a/tests/li_solv/reference/li.json b/tests/li_solv/reference/li.json
index 2a4a75618..fde26d561 100644
--- a/tests/li_solv/reference/li.json
+++ b/tests/li_solv/reference/li.json
@@ -14,15 +14,17 @@
     },
     "molecule": {
       "cavity": {
-        "sigma": 0.5,
         "spheres": [
           {
+            "alpha": 1.0,
+            "beta": 0.0,
             "center": [
               0.0,
               0.0,
               0.0
             ],
-            "radius": 4.0
+            "radius": 4.0,
+            "sigma": 0.5
           }
         ]
       },
@@ -30,6 +32,7 @@
       "coords": [
         {
           "atom": "li",
+          "r_rms": 4.0992133976e-05,
           "xyz": [
             0.0,
             0.0,
@@ -83,19 +86,19 @@
           "derivative": "abgv_55"
         },
         "nuclear_operator": {
+          "nuclear_model": "point_like",
           "proj_prec": 0.001,
           "shared_memory": false,
           "smooth_prec": 0.001
         },
         "reaction_operator": {
           "density_type": "total",
-          "dynamic_thrs": true,
+          "dynamic_thrs": false,
           "epsilon_in": 1.0,
           "epsilon_out": 2.0,
           "formulation": "exponential",
           "kain": 5,
           "max_iter": 100,
-          "optimizer": "potential",
           "poisson_prec": 0.001
         },
         "xc_operator": {
@@ -162,7 +165,7 @@
       "charge": 1,
       "dipole_moment": {
         "dip-1": {
-          "magnitude": 7.201515824859437e-14,
+          "magnitude": 7.400866885254787e-14,
           "r_O": [
             0.0,
             0.0,
@@ -198,7 +201,7 @@
       "multiplicity": 1,
       "orbital_energies": {
         "energy": [
-          -2.2040278701953486
+          -2.204016950767741
         ],
         "occupation": [
           2.0
@@ -206,23 +209,23 @@
         "spin": [
           "p"
         ],
-        "sum_occupied": -4.408055740390697
+        "sum_occupied": -4.408033901535482
       },
       "scf_energy": {
-        "E_ee": 3.332799383894309,
+        "E_ee": 3.3327993838943106,
         "E_eext": 0.0,
-        "E_el": -7.094583823998505,
-        "E_en": -16.53075528769529,
-        "E_kin": 7.675564301900199,
+        "E_el": -7.094572904570907,
+        "E_en": -16.530755287695293,
+        "E_kin": 7.675564301900204,
         "E_next": 0.0,
         "E_nn": 0.0,
-        "E_nuc": -0.1921541203648782,
-        "E_tot": -7.2867379443633835,
-        "E_x": -0.41648449531759685,
-        "E_xc": -1.2837729977609404,
-        "Er_el": 0.12806527098081377,
-        "Er_nuc": -0.1921541203648782,
-        "Er_tot": -0.06408884938406478
+        "E_nuc": -0.19217053480332283,
+        "E_tot": -7.28674343937423,
+        "E_x": -0.4164844953175972,
+        "E_xc": -1.2837729977609416,
+        "Er_el": 0.12807619040840895,
+        "Er_nuc": -0.19217053480332283,
+        "Er_tot": -0.06409434439491388
       }
     },
     "provenance": {
@@ -231,25 +234,25 @@
       "nthreads": 1,
       "routine": "mrchem.x",
       "total_cores": 1,
-      "version": "1.1.0-alpha"
+      "version": "1.2.0-alpha"
     },
     "rsp_calculations": null,
     "scf_calculation": {
       "initial_energy": {
-        "E_ee": 3.332799383894309,
+        "E_ee": 3.3327993838943106,
         "E_eext": 0.0,
-        "E_el": -7.094583823998505,
-        "E_en": -16.53075528769529,
-        "E_kin": 7.675564301900199,
+        "E_el": -7.094572904570907,
+        "E_en": -16.530755287695293,
+        "E_kin": 7.675564301900204,
         "E_next": 0.0,
         "E_nn": 0.0,
-        "E_nuc": -0.1921541203648782,
-        "E_tot": -7.2867379443633835,
-        "E_x": -0.41648449531759685,
-        "E_xc": -1.2837729977609404,
-        "Er_el": 0.12806527098081377,
-        "Er_nuc": -0.1921541203648782,
-        "Er_tot": -0.06408884938406478
+        "E_nuc": -0.19217053480332283,
+        "E_tot": -7.28674343937423,
+        "E_x": -0.4164844953175972,
+        "E_xc": -1.2837729977609416,
+        "Er_el": 0.12807619040840895,
+        "Er_nuc": -0.19217053480332283,
+        "Er_tot": -0.06409434439491388
       },
       "success": true
     },
diff --git a/tests/li_solv/reference/li.out b/tests/li_solv/reference/li.out
index 6512997f6..4f21e9988 100644
--- a/tests/li_solv/reference/li.out
+++ b/tests/li_solv/reference/li.out
@@ -9,12 +9,12 @@
 ***         | |  | |  _ <| |___| | | |  __/ | | | | |                   ***
 ***         |_|  |_|_| \_\\____|_| |_|\___|_| |_| |_|                   ***
 ***                                                                     ***
-***         VERSION            1.1.0-alpha                              ***
+***         VERSION            1.2.0-alpha                              ***
 ***                                                                     ***
-***         Git branch         solvent-input                            ***
-***         Git commit hash    685b3596305e75a76eff-dirty               ***
-***         Git commit author  Stig Rune Jensen                         ***
-***         Git commit date    Wed Jul 13 16:03:41 2022 +0200           ***
+***         Git branch         vecmult                                  ***
+***         Git commit hash    07cbcd62711b3df44973-dirty               ***
+***         Git commit author  Gabriel Gerez                            ***
+***         Git commit date    Mon Oct 30 13:15:37 2023 +0100           ***
 ***                                                                     ***
 ***         Contact: luca.frediani@uit.no                               ***
 ***                                                                     ***
@@ -46,11 +46,11 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
 
 ---------------------------------------------------------------------------
 
- MRCPP version         : 1.4.1
+ MRCPP version         : 1.6.0-alpha
  Git branch            : HEAD
- Git commit hash       : 75d41879b1908a94a452
- Git commit author     : Stig Rune Jensen
- Git commit date       : Thu Jan 6 11:38:53 2022 +0100
+ Git commit hash       : f8def0a086da6410e5dd
+ Git commit author     : Peter Wind
+ Git commit date       : Sat Oct 21 13:32:25 2023 +0200
 
  Linear algebra        : EIGEN v3.4.0
  Parallelization       : OpenMP (1 threads)
@@ -101,13 +101,12 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
                              Solvation Cavity
 ---------------------------------------------------------------------------
  Formulation             :                                     exponential
- Cavity width            :                                        0.500000
  Dielectric constant     :            (in)                        1.000000
                          :           (out)                        2.000000
 ---------------------------------------------------------------------------
-    N      Radius        :               x               y               z
+    N      R_0  Alpha Beta Sigma      Radius         x         y         z
 ---------------------------------------------------------------------------
-    0    4.000000        :        0.000000        0.000000        0.000000
+    0   4.0000  1.00  0.00  0.50  ->  4.0000  0.000000  0.000000  0.000000
 ===========================================================================
 
 
@@ -165,17 +164,17 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
  X-C energy              :            (au)                 -1.283772997761
  N-N energy              :            (au)                  0.000000000000
 ---------------------------------------------------------------------------
- Reaction energy (el)    :            (au)                  0.128065270981
- Reaction energy (nuc)   :            (au)                 -0.192154120365
- Reaction energy (tot)   :            (au)                 -0.064088849384
+ Reaction energy (el)    :            (au)                  0.128076190408
+ Reaction energy (nuc)   :            (au)                 -0.192170534803
+ Reaction energy (tot)   :            (au)                 -0.064094344395
 ---------------------------------------------------------------------------
- Electronic energy       :            (au)                 -7.094583823999
- Nuclear energy          :            (au)                 -0.192154120365
+ Electronic energy       :            (au)                 -7.094572904571
+ Nuclear energy          :            (au)                 -0.192170534803
 ---------------------------------------------------------------------------
- Total energy            :            (au)             -7.286737944363e+00
-                         :      (kcal/mol)             -4.572497095103e+03
-                         :        (kJ/mol)             -1.913132784591e+04
-                         :            (eV)             -1.982822406774e+02
+ Total energy            :            (au)             -7.286743439374e+00
+                         :      (kcal/mol)             -4.572500543274e+03
+                         :        (kJ/mol)             -1.913134227306e+04
+                         :            (eV)             -1.982823902042e+02
 ===========================================================================
 
 
@@ -184,9 +183,9 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
 ---------------------------------------------------------------------------
     n  Occ Spin          :                                         Epsilon
 ---------------------------------------------------------------------------
-    0    2    p          :            (au)                 -2.204027870195
+    0    2    p          :            (au)                 -2.204016950768
 ---------------------------------------------------------------------------
- Sum occupied            :            (au)                 -4.408055740391
+ Sum occupied            :            (au)                 -4.408033901535
 ===========================================================================
 
 
@@ -222,17 +221,17 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
  X-C energy              :            (au)                 -1.283772997761
  N-N energy              :            (au)                  0.000000000000
 ---------------------------------------------------------------------------
- Reaction energy (el)    :            (au)                  0.128065270981
- Reaction energy (nuc)   :            (au)                 -0.192154120365
- Reaction energy (tot)   :            (au)                 -0.064088849384
+ Reaction energy (el)    :            (au)                  0.128076190408
+ Reaction energy (nuc)   :            (au)                 -0.192170534803
+ Reaction energy (tot)   :            (au)                 -0.064094344395
 ---------------------------------------------------------------------------
- Electronic energy       :            (au)                 -7.094583823999
- Nuclear energy          :            (au)                 -0.192154120365
+ Electronic energy       :            (au)                 -7.094572904571
+ Nuclear energy          :            (au)                 -0.192170534803
 ---------------------------------------------------------------------------
- Total energy            :            (au)             -7.286737944363e+00
-                         :      (kcal/mol)             -4.572497095103e+03
-                         :        (kJ/mol)             -1.913132784591e+04
-                         :            (eV)             -1.982822406774e+02
+ Total energy            :            (au)             -7.286743439374e+00
+                         :      (kcal/mol)             -4.572500543274e+03
+                         :        (kJ/mol)             -1.913134227306e+04
+                         :            (eV)             -1.982823902042e+02
 ===========================================================================
 
 
@@ -241,9 +240,9 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
 ---------------------------------------------------------------------------
     n  Occ Spin          :                                         Epsilon
 ---------------------------------------------------------------------------
-    0    2    p          :            (au)                 -2.204027870195
+    0    2    p          :            (au)                 -2.204016950768
 ---------------------------------------------------------------------------
- Sum occupied            :            (au)                 -4.408055740391
+ Sum occupied            :            (au)                 -4.408033901535
 ===========================================================================
 
 
@@ -272,7 +271,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
 ***                                                                     ***
 ***                            Exiting MRChem                           ***
 ***                                                                     ***
-***                       Wall time :  0h  0m 51s                       ***
+***                       Wall time :  0h  0m 53s                       ***
 ***                                                                     ***
 ***************************************************************************
                                                                            
diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp
index 9468e465b..045494098 100644
--- a/tests/solventeffect/reaction_operator.cpp
+++ b/tests/solventeffect/reaction_operator.cpp
@@ -38,10 +38,12 @@
 #include "chemistry/Element.h"
 #include "chemistry/Nucleus.h"
 #include "chemistry/PeriodicTable.h"
+#include "chemistry/chemistry_utils.h"
 #include "environment/Cavity.h"
 #include "environment/Permittivity.h"
 #include "environment/SCRF.h"
 #include "qmfunctions/Orbital.h"
+#include "qmfunctions/density_utils.h"
 #include "qmfunctions/orbital_utils.h"
 #include "qmoperators/two_electron/ReactionOperator.h"
 
@@ -85,14 +87,22 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") {
     HydrogenFunction f(1, 0, 0);
     if (mrcpp::mpi::my_orb(Phi[0])) mrcpp::cplxfunc::project(Phi[0], f, NUMBER::Real, prec);
 
+    auto rho_nuc = chemistry::compute_nuclear_density(prec, molecule, 100);
+
     int kain = 4;
-    auto scrf_p = std::make_unique<SCRF>(dielectric_func, molecule, P_p, D_p, prec, kain, 100, true, false, "total");
+    auto scrf_p = std::make_unique<SCRF>(dielectric_func, rho_nuc, P_p, D_p, kain, 100, false, "total");
 
     auto Reo = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
     Reo->setup(prec);
-    double total_energy = Reo->getTotalEnergy();
+
+    Density rho_el(false);
+    density::compute(prec, rho_el, Phi, DensityType::Total);
+    rho_el.rescale(-1.0);
+
+    auto [Er_nuc, Er_el] = Reo->getHelper()->computeEnergies(rho_el);
+    auto total_energy = Er_nuc + Er_el;
     Reo->clear();
-    REQUIRE(total_energy == Approx(-0.191434124263).epsilon(thrs));
+    REQUIRE(total_energy == Approx(-1.022729683846e-01).epsilon(thrs));
 }
 
 } // namespace reaction_operator