diff --git a/src/environment/GPESolver.cpp b/src/environment/GPESolver.cpp
index df7141c04..f27937df7 100644
--- a/src/environment/GPESolver.cpp
+++ b/src/environment/GPESolver.cpp
@@ -121,17 +121,12 @@ void GPESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFu
 mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, const Density &rho_el) {
     mrcpp::ComplexFunction Poisson_func;
     mrcpp::ComplexFunction rho_eff;
-    mrcpp::ComplexFunction first_term;
     mrcpp::ComplexFunction Vr_np1;
     Vr_np1.alloc(NUMBER::Real);
 
-    auto eps_inv_func = mrcpp::AnalyticFunction<3>([this](const mrcpp::Coord<3> &r) { return 1.0 / this->epsilon.evalf(r); });
     Density rho_tot(false);
     computeDensities(rho_el, rho_tot);
-
-    mrcpp::cplxfunc::multiply(first_term, rho_tot, eps_inv_func, this->apply_prec);
-
-    mrcpp::cplxfunc::add(rho_eff, 1.0, first_term, -1.0, rho_tot, -1.0);
+    computeEffectiveDensity(rho_eff, rho_tot);
     rho_tot.free(NUMBER::Real);
 
     mrcpp::cplxfunc::add(Poisson_func, 1.0, in_gamma, 1.0, rho_eff, -1.0);
@@ -139,6 +134,16 @@ mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunct
     return Vr_np1;
 }
 
+void GPESolver::computeEffectiveDensity(mrcpp::ComplexFunction &rho_eff, Density &rho_tot) const {
+    mrcpp::ComplexFunction rho_eps;
+
+    auto eps_inv_func = mrcpp::AnalyticFunction<3>([this](const mrcpp::Coord<3> &r) { return 1.0 / this->epsilon.evalf(r); });
+
+    mrcpp::cplxfunc::multiply(rho_eps, rho_tot, eps_inv_func, this->apply_prec);
+
+    mrcpp::cplxfunc::add(rho_eff, 1.0, rho_eps, -1.0, rho_tot, -1.0);
+}
+
 void GPESolver::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain) {
     OrbitalVector phi_n(0);
     OrbitalVector dPhi_n(0);
@@ -190,6 +195,10 @@ void GPESolver::runMicroIterations(const mrcpp::ComplexFunction &V_vac, const De
             Vr_np1.free(NUMBER::Real);
             mrcpp::cplxfunc::add(Vr_np1, 1.0, Vr_n, 1.0, dVr_n, -1.0);
         }
+        // DEBUG
+        auto [Er_el, Er_nuc] = computeEnergies(rho_el);
+        std::cout << "\niter: " << iter << "\nenergy: " << Er_nuc << "\nupdate: " << update << std::endl;
+        // DEBUG
 
         // set up for next iteration
         resetComplexFunction(this->Vr_n);
diff --git a/src/environment/GPESolver.h b/src/environment/GPESolver.h
index acd8adf8e..181f0b98a 100644
--- a/src/environment/GPESolver.h
+++ b/src/environment/GPESolver.h
@@ -161,6 +161,8 @@ class GPESolver {
      */
     virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma);
 
+    virtual void computeEffectiveDensity(mrcpp::ComplexFunction &rho_eff, Density &rho_tot) const;
+
     /** @brief Iterates once through the Generalized Poisson equation to compute the reaction potential
      * @param ingamma the surface charge distribution
      * @param Phi the molecular orbitals
diff --git a/src/environment/PBESolver.cpp b/src/environment/PBESolver.cpp
index 6c5e22fd6..f6990a27b 100644
--- a/src/environment/PBESolver.cpp
+++ b/src/environment/PBESolver.cpp
@@ -96,4 +96,21 @@ void PBESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFu
     out_gamma.add(-1.0, pb_term);
 }
 
+void PBESolver::computeEffectiveDensity(mrcpp::ComplexFunction &rho_eff, Density &rho_tot) const {
+    mrcpp::ComplexFunction rho_eps;
+
+    auto eps_inv_func = mrcpp::AnalyticFunction<3>([this](const mrcpp::Coord<3> &r) { return 1.0 / this->epsilon.evalf(r); });
+
+    const double k_b = 1.380649e-23;
+    const double T = 298.15;
+    const double J2H = 1 / (4.3597447222071e-18);
+    double factor = 1 / (k_b * T * J2H);
+
+    std::cout << "factor: " << factor << std::endl;
+
+    mrcpp::cplxfunc::multiply(rho_eps, rho_tot, eps_inv_func, this->apply_prec);
+
+    mrcpp::cplxfunc::add(rho_eff, factor, rho_eps, -1.0, rho_tot, -1.0);
+}
+
 } // namespace mrchem
diff --git a/src/environment/PBESolver.h b/src/environment/PBESolver.h
index c5215d55b..e3c10bd80 100644
--- a/src/environment/PBESolver.h
+++ b/src/environment/PBESolver.h
@@ -90,5 +90,7 @@ class PBESolver : public GPESolver {
      * @details The PB term is computed as \f$ \kappa^2 \sinh(V_{tot}) \f$ and returned.
      */
     virtual void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term);
+
+    void computeEffectiveDensity(mrcpp::ComplexFunction &rho_eff, Density &rho_tot) const override;
 };
 } // namespace mrchem
diff --git a/tests/solventeffect/PB_solver.cpp b/tests/solventeffect/PB_solver.cpp
index 1721b6662..35faa26ed 100644
--- a/tests/solventeffect/PB_solver.cpp
+++ b/tests/solventeffect/PB_solver.cpp
@@ -53,7 +53,7 @@ namespace PB_solver {
  */
 
 TEST_CASE("Poisson Boltzmann equation solver standard", "[PB_solver][pb_standard]") {
-    const double prec = 1.0e-3;
+    const double prec = 1.0e-5;
     const double thrs = 1.0e-4;
 
     auto dyn_thrs = false;
@@ -62,7 +62,7 @@ TEST_CASE("Poisson Boltzmann equation solver standard", "[PB_solver][pb_standard
     auto eps_in = 1.0;
     auto eps_out = 78.54;
     auto kappa_out = 0.054995;
-    auto slope = 0.1;
+    auto slope = 0.2;
 
     auto R = std::vector<double>({3.7794522509156563});
     auto sph_coords = std::vector<mrcpp::Coord<3>>({{0.0, 0.0, 0.0}});
@@ -103,12 +103,13 @@ TEST_CASE("Poisson Boltzmann equation solver standard", "[PB_solver][pb_standard
 
         auto [Er_el, Er_nuc] = Reo->getSolver()->computeEnergies(rho_el);
         Reo->clear();
+        std::cout << "Er_nuc: " << Er_nuc << std::endl;
         REQUIRE((Er_nuc) == Approx(-1.329978908155e-01).epsilon(thrs)); // exact is -0.1373074208 Hartree, though ours is close, i think we are a bit too far away, some parameterization issue
     }
 }
 
 TEST_CASE("Poisson Boltzmann equation solver linearized", "[PB_solver][pb_linearized]") {
-    const double prec = 1.0e-3;
+    const double prec = 1.0e-5;
     const double thrs = 1.0e-4;
 
     auto dyn_thrs = false;
@@ -117,7 +118,7 @@ TEST_CASE("Poisson Boltzmann equation solver linearized", "[PB_solver][pb_linear
     auto eps_in = 1.0;
     auto eps_out = 78.54;
     auto kappa_out = 0.054995;
-    auto slope = 0.1;
+    auto slope = 0.2;
 
     auto R = std::vector<double>({3.7794522509156563});
     auto sph_coords = std::vector<mrcpp::Coord<3>>({{0.0, 0.0, 0.0}});
@@ -160,7 +161,10 @@ TEST_CASE("Poisson Boltzmann equation solver linearized", "[PB_solver][pb_linear
 
         auto [Er_el, Er_nuc] = Reo->getSolver()->computeEnergies(rho_el);
         Reo->clear();
+        std::cout << "Er_nuc: " << Er_nuc << std::endl;
         REQUIRE(Er_nuc == Approx(-1.329978908155e-01).epsilon(thrs)); // what we get in standard GPESolver is -1.455145361712e-01, while with PB we get -1.329978908155e-01
+        // Er_nuc: -1.585665435733e-01
+        // 50: Er_nuc: -1.585671452611e-01
     }
 }