From 00af3001587a1e85e21443221a65492b2b7dc723 Mon Sep 17 00:00:00 2001 From: Gabriel Gerez Date: Tue, 26 Sep 2023 16:10:26 +0200 Subject: [PATCH] Rename GPESolver::solveEquation to iterateEquation Remove print statements Make sure the code compiles --- pilot/mrchem.cpp | 2 +- src/driver.cpp | 1 - src/environment/GPESolver.cpp | 42 ++++--------------- src/environment/GPESolver.h | 2 +- src/environment/LPBESolver.h | 2 - src/environment/PBESolver.h | 6 ++- .../two_electron/ReactionPotential.cpp | 2 +- tests/solventeffect/reaction_operator.cpp | 1 + 8 files changed, 15 insertions(+), 43 deletions(-) diff --git a/pilot/mrchem.cpp b/pilot/mrchem.cpp index 063f9b71b..fe449ec97 100644 --- a/pilot/mrchem.cpp +++ b/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(); diff --git a/src/driver.cpp b/src/driver.cpp index 56f53bcfa..112402972 100644 --- a/src/driver.cpp +++ b/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(); diff --git a/src/environment/GPESolver.cpp b/src/environment/GPESolver.cpp index 95e73cbbd..3cdd2966a 100644 --- a/src/environment/GPESolver.cpp +++ b/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); diff --git a/src/environment/GPESolver.h b/src/environment/GPESolver.h index 6a567cbd7..59e2767cb 100644 --- a/src/environment/GPESolver.h +++ b/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 &Phi); + mrcpp::ComplexFunction &iterateEquation(double prec, const std::shared_ptr &Phi); double getTotalEnergy(); void clear(); diff --git a/src/environment/LPBESolver.h b/src/environment/LPBESolver.h index 7c0ef21ee..59edb8a2e 100644 --- a/src/environment/LPBESolver.h +++ b/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: diff --git a/src/environment/PBESolver.h b/src/environment/PBESolver.h index 4ae37dbd3..4489af802 100644 --- a/src/environment/PBESolver.h +++ b/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); }; diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp index 24924f980..331bbfd62 100644 --- a/src/qmoperators/two_electron/ReactionPotential.cpp +++ b/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); diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp index d96120d52..5c86c55d4 100644 --- a/tests/solventeffect/reaction_operator.cpp +++ b/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)); }