Rename GPESolver::solveEquation to iterateEquation

Remove print statements
Make sure the code compiles

pilot/mrchem.cpp

@@ -125,7 +125,7 @@ int main(int argc, char **argv) {
 
     PB_solver->setConvergenceThreshold(poisson_prec);
     // PB_solver->setStaticSalt();
-    PB_solver->solveEquation(poisson_prec, ORB);
+    PB_solver->iterateEquation(poisson_prec, ORB);
     double total_energy = PB_solver->getTotalEnergy() / 2;
     PB_solver->clear();
 

src/driver.cpp

@@ -1066,7 +1066,6 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
 
         auto width_ion = json_fock["reaction_operator"]["ion_width"];
         auto kformulation = json_fock["reaction_operator"]["formulation"];
-        auto solver_type = json_fock["reaction_operator"]["solver_type"];
 
         Permittivity dielectric_func(*cavity_p, eps_i, eps_o, formulation);
         dielectric_func.printParameters();

src/environment/GPESolver.cpp

@@ -71,8 +71,6 @@ GPESolver::GPESolver(Permittivity e, const Nuclei &N, PoissonOperator_p P, Deriv
         , poisson(P) {
     setDCavity();
     rho_nuc = chemistry::compute_nuclear_density(this->apply_prec, N, 1000);
-    std::cout << "norm of rho_nuc: " << rho_nuc.norm() << std::endl;
-    std::cout << "integral of rho_nuc: " << rho_nuc.integrate() << std::endl;
 }
 
 GPESolver::~GPESolver() {
@@ -115,10 +113,10 @@ void GPESolver::computeDensities(OrbitalVector &Phi) {
     rho_el.rescale(-1.0);
     if (this->density_type == "electronic") {
         mrcpp::cplxfunc::deep_copy(this->rho_tot, rho_el);
+
     } else if (this->density_type == "nuclear") {
         mrcpp::cplxfunc::deep_copy(this->rho_tot, this->rho_nuc);
-        std::cout << "norm of rho_tot: " << this->rho_tot.norm() << std::endl;
-        std::cout << "integral of rho_tot: " << this->rho_tot.integrate() << std::endl;
+
     } else {
         mrcpp::cplxfunc::add(this->rho_tot, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); // probably change this into a vector
     }
@@ -131,8 +129,6 @@ void GPESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFu
     mrcpp::dot(this->apply_prec, out_gamma.real(), d_V, this->d_cavity);
     out_gamma.rescale(std::log((epsilon.getEpsIn() / epsilon.getEpsOut())) * (1.0 / (4.0 * mrcpp::pi)));
     mrcpp::clear(d_V, true);
-    std::cout << "norm of gamma: " << out_gamma.norm() << std::endl;
-    std::cout << "integral of gamma: " << out_gamma.integrate() << std::endl;
 }
 
 mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma) {
@@ -150,9 +146,6 @@ mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunct
     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);
 
-    std::cout << "norm of rho_eff: " << rho_eff.norm() << std::endl;
-    std::cout << "integral of rho_eff: " << rho_eff.integrate() << std::endl;
-
     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());
 
@@ -185,15 +178,12 @@ void GPESolver::runMicroIterations(mrcpp::ComplexFunction V_vac) {
     auto update = 10.0, norm = 1.0;
 
     auto iter = 1;
-    std::cout << "convergence threshold: " << this->conv_thrs << std::endl;
     for (; iter <= max_iter; iter++) {
 
         Timer t_iter;
         mrcpp::ComplexFunction Vr_np1;
         // solve the poisson equation
         Vr_np1 = solvePoissonEquation(this->gamma_n);
-        std::cout << "norm of Vr_np1: " << Vr_np1.norm() << std::endl;
-        std::cout << "integral of Vr_np1: " << Vr_np1.integrate() << std::endl;
 
         norm = Vr_np1.norm();
 
@@ -207,7 +197,6 @@ void GPESolver::runMicroIterations(mrcpp::ComplexFunction V_vac) {
             mrcpp::cplxfunc::add(Vr_np1, 1.0, Vr_n, 1.0, dVr_n, -1.0);
         }
         update = dVr_n.norm();
-        std::cout << "total energy at iteration " << iter << ": " << getTotalEnergy() << std::endl;
         // set up for next iteration
         mrcpp::ComplexFunction V_tot;
         mrcpp::cplxfunc::add(V_tot, 1.0, Vr_np1, 1.0, V_vac, -1.0);
@@ -229,9 +218,7 @@ void GPESolver::runMicroIterations(mrcpp::ComplexFunction V_vac) {
         updateCurrentGamma(gamma_np1);
         printConvergenceRow(iter, norm, update, t_iter.elapsed());
 
-        std::cout << "update: " << update << std::endl;
         // compute new PB term
-        std::cout << "update " << update << " convergence threshold " << this->conv_thrs << std::endl;
         if (update < this->conv_thrs) break;
     }
 
@@ -268,15 +255,13 @@ void GPESolver::printConvergenceRow(int i, double norm, double update, double ti
     println(3, o_txt.str());
 }
 
-mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const OrbitalVector_p &Phi) {
+mrcpp::ComplexFunction &GPESolver::iterateEquation(double prec, const OrbitalVector_p &Phi) {
     this->apply_prec = prec;
     computeDensities(*Phi);
     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());
-    std::cout << "norm of V_vac: " << V_vac.norm() << std::endl;
-    std::cout << "integral of V_vac: " << V_vac.integrate() << std::endl;
     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
@@ -286,17 +271,13 @@ mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const OrbitalVecto
         mrcpp::ComplexFunction gamma_0;
         mrcpp::ComplexFunction V_tot;
         computeGamma(V_vac, gamma_0);
-        std::cout << "norm of gamma_n at zeroth iteration: " << gamma_0.norm() << std::endl;
-        std::cout << "integral of gamma_n at zeroth iteration: " << gamma_0.integrate() << std::endl;
+
         this->Vr_n = solvePoissonEquation(gamma_0);
-        std::cout << "total energy at zeroth iteration" << getTotalEnergy() << std::endl;
-        std::cout << "norm of Vr_n at zeroth iteration: " << this->Vr_n.norm() << std::endl;
-        std::cout << "integral of Vr_n at zeroth iteration: " << this->Vr_n.integrate() << std::endl;
+
         mrcpp::cplxfunc::add(V_tot, 1.0, V_vac, 1.0, this->Vr_n, -1.0);
         computeGamma(V_tot, this->gamma_n);
     }
-    std::cout << "norm of gamma_n at first iteration: " << this->gamma_n.norm() << std::endl;
-    std::cout << "integral of gamma_n at first iteration: " << this->gamma_n.integrate() << std::endl;
+
     // update the potential/gamma with the previous scf cycle kain update before doing anything with them
 
     if (accelerate_Vr) {
@@ -308,16 +289,7 @@ mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const OrbitalVecto
         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);
-        // std::cout << "apply prec " << this->apply_prec << std::endl;
-        // std::cout << "conv thrs " << this->conv_thrs << std::endl;
-        // std::cout << "difference between them " << this->conv_thrs - this->apply_prec << std::endl;
-        // if (std::abs(this->conv_thrs - this->apply_prec) <= mrcpp::MachineZero) {
-        //     std::cout  << "does this ever work?" << std::endl;
-        //     this->salt_factor = 1.0;
-        //     computePBTerm(V_tot, this->salt_factor);
-        // } else {
-        //     this->salt_factor = 0.0;
-        // }
+
     } else {
         mrcpp::ComplexFunction temp_gamma_n;
         mrcpp::cplxfunc::add(temp_gamma_n, 1.0, this->gamma_n, 1.0, this->dgamma_n, -1.0);

src/environment/GPESolver.h

@@ -80,7 +80,7 @@ class GPESolver {
     void updateMOResidual(double const err_t) { this->mo_residual = err_t; }
 
     friend class ReactionPotential;
-    mrcpp::ComplexFunction &solveEquation(double prec, const std::shared_ptr<mrchem::OrbitalVector> &Phi);
+    mrcpp::ComplexFunction &iterateEquation(double prec, const std::shared_ptr<mrchem::OrbitalVector> &Phi);
     double getTotalEnergy();
 
     void clear();

src/environment/LPBESolver.h

@@ -52,8 +52,6 @@ class LPBESolver final : public PBESolver {
                std::string density_type)
             : PBESolver(e, kappa, N, P, D, orb_prec, kain_hist, max_iter, acc_pot, dyn_thrs, density_type) {}
 
-    ~LPBESolver();
-
     friend class ReactionPotential;
 
 protected:

src/environment/PBESolver.h

@@ -51,8 +51,6 @@ class PBESolver : public GPESolver {
               bool dyn_thrs,
               std::string density_type);
 
-    ~PBESolver();
-
     void setStaticSalt(bool do_static) { do_static_salt = do_static; }
 
     friend class ReactionPotential;
@@ -65,6 +63,10 @@ class PBESolver : public GPESolver {
 
     void setKappa(DHScreening k);
 
+    // FIXME ComputeGamma should not be computing the PB term.
+    //  THe PB term is actually an approximation for an external density, so
+    //  it should be computed in the computedensities function or in the
+    //  solvePoissonEquation function.
     void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma);
     void computePBTerm(mrcpp::ComplexFunction &V_tot, double salt_factor);
 };

src/qmoperators/two_electron/ReactionPotential.cpp

@@ -49,7 +49,7 @@ void ReactionPotential::setup(double prec) {
     mrcpp::print::value(3, "Precision", prec, "(rel)", 5);
     mrcpp::print::value(3, "Threshold", thrs, "(abs)", 5);
     mrcpp::print::separator(3, '-');
-    auto potential = this->solver->solveEquation(prec, this->Phi);
+    auto potential = this->solver->iterateEquation(prec, this->Phi);
     mrcpp::cplxfunc::deep_copy(*this, potential);
     if (plevel == 2) print_utils::qmfunction(2, "Reaction operator", *this, timer);
     mrcpp::print::footer(3, timer, 2);

tests/solventeffect/reaction_operator.cpp

@@ -92,6 +92,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") {
     Reo->setup(prec);
     double total_energy = Reo->getTotalEnergy();
     Reo->clear();
+    std::cout << "total_energy = " << total_energy << std::endl;
     REQUIRE(total_energy == Approx(-0.191434124263).epsilon(thrs));
 }