fix tests for pb solver

MRChemSoft · Nov 14, 2023 · 5a3f77d · 5a3f77d
1 parent b20a718
commit 5a3f77d
Show file tree

Hide file tree

Showing 21 changed files with 12,375 additions and 19,619 deletions.
diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst
diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py
diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst
diff --git a/src/driver.cpp b/src/driver.cpp
@@ -1070,7 +1070,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         dielectric_func.printParameters();
         Density rho_nuc(false);
         rho_nuc = chemistry::compute_nuclear_density(poisson_prec, nuclei, 100);
-        std::shared_ptr<GPESolver> scrf_p;
+        std::unique_ptr<GPESolver> scrf_p;
         auto solver_type = json_fock["reaction_operator"]["solver_type"];
 
         if (solver_type == "Poisson-Boltzmann") {
@@ -1097,7 +1097,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
             dhscreening.printParameters();
             scrf_p = std::make_unique<LPBESolver>(dielectric_func, dhscreening, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type);
         } else if (solver_type == "Generalized_Poisson") {
-            scrf_p = auto scrf_p = std::make_unique<GPESolver>(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type);
+            scrf_p = std::make_unique<GPESolver>(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type);
         } else {
             MSG_ERROR("Solver type not implemented");
         }

diff --git a/src/environment/GPESolver.cpp b/src/environment/GPESolver.cpp
@@ -195,6 +195,7 @@ void GPESolver::runMicroIterations(const mrcpp::ComplexFunction &V_vac, std::sha
 
         // compute new PB term
         if (update < this->conv_thrs) break;
+        iter++;
     }
 
     if (iter > max_iter) println(0, "Reaction potential failed to converge after " << iter - 1 << " iterations, residual " << update);
@@ -285,7 +286,7 @@ void GPESolver::printParameters() const {
     }
 
     nlohmann::json data = {
-        {"Method                ", "GPE Solver"},
+        {"Method                ", this->solver_name},
         {"Density               ", this->density_type},
         {"Max iterations        ", max_iter},
         {"KAIN solver           ", o_kain.str()},

diff --git a/src/environment/GPESolver.h b/src/environment/GPESolver.h
@@ -59,13 +59,12 @@ class GPESolver {
 
     friend class ReactionPotential;
     auto computeEnergies(const Density &rho_el) -> std::tuple<double, double>;
-    void clear();
 
 protected:
     void clear();
-
     bool dynamic_thrs;
     std::string density_type;
+    std::string solver_name = "Generalized Poisson";
 
     int max_iter;
     int history;
@@ -100,7 +99,7 @@ class GPESolver {
 
     void resetComplexFunction(mrcpp::ComplexFunction &function);
 
-    void printParameters() const;
+    virtual void printParameters() const;
     void printConvergenceRow(int i, double norm, double update, double time) const;
 };
 } // namespace mrchem
diff --git a/src/environment/LPBESolver.cpp b/src/environment/LPBESolver.cpp
@@ -47,12 +47,22 @@ using DerivativeOperator_p = std::shared_ptr<mrcpp::DerivativeOperator<3>>;
 using OrbitalVector_p = std::shared_ptr<mrchem::OrbitalVector>;
 
 namespace mrchem {
-
+LPBESolver::LPBESolver(const Permittivity &e,
+                       const DHScreening &k,
+                       const Density &rho_nuc,
+                       std::shared_ptr<mrcpp::PoissonOperator> P,
+                       std::shared_ptr<mrcpp::DerivativeOperator<3>> D,
+                       int kain_hist,
+                       int max_iter,
+                       bool dyn_thrs,
+                       const std::string &density_type)
+        : PBESolver(e, k, rho_nuc, P, D, kain_hist, max_iter, dyn_thrs, density_type)
+        , solver_name("Linearized Poisson-Boltzmann") {}
 // TODO separate this for the linear and non-linear solver
-void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) {
+void LPBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) {
     resetComplexFunction(pb_term);
-    mrcpp::cplxfunc::multiply(pb_term, this->kappa, V_tot, this->apply_prec);
-    pbe_term.rescale(salt_factor);
+    mrcpp::cplxfunc::multiply(pb_term, V_tot, this->kappa, this->apply_prec);
+    pb_term.rescale(salt_factor / (4.0 * mrcpp::pi));
 }
 
 } // namespace mrchem
diff --git a/src/environment/LPBESolver.h b/src/environment/LPBESolver.h
@@ -47,12 +47,12 @@ class LPBESolver final : public PBESolver {
                int kain_hist,
                int max_iter,
                bool dyn_thrs,
-               const std::string &density_type)
-            : PBESolver(e, kappa, N, P, D, orb_prec, kain_hist, max_iter, acc_pot, dyn_thrs, density_type) {}
+               const std::string &density_type);
 
     friend class ReactionPotential;
 
 protected:
     void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) override;
+    std::string solver_name;
 };
 } // namespace mrchem
diff --git a/src/environment/PBESolver.cpp b/src/environment/PBESolver.cpp
@@ -58,7 +58,8 @@ PBESolver::PBESolver(const Permittivity &e,
                      bool dyn_thrs,
                      const std::string &density_type)
         : GPESolver(e, rho_nuc, P, D, kain_hist, max_iter, dyn_thrs, density_type)
-        , kappa(k) {}
+        , kappa(k)
+        , solver_name("Poisson-Boltzmann") {}
 
 void PBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) {
     // create a lambda function for the sinh(V) term and multiply it with kappa and salt factor to get the PB term
@@ -68,7 +69,7 @@ void PBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_f
     sinhV.alloc(NUMBER::Real);
     mrcpp::map(this->apply_prec / 100, sinhV.real(), V_tot.real(), sinh_f);
 
-    mrcpp::cplxfunc::multiply(pb_term, this->kappa, sinhV, this->apply_prec);
+    mrcpp::cplxfunc::multiply(pb_term, sinhV, this->kappa, this->apply_prec);
 }
 
 void PBESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) {

diff --git a/src/environment/PBESolver.h b/src/environment/PBESolver.h
@@ -54,6 +54,7 @@ class PBESolver : public GPESolver {
 protected:
     DHScreening kappa;
     mrcpp::ComplexFunction rho_ext;
+    std::string solver_name;
 
     // FIXME ComputeGamma should not be computing the PB term.
     //  THe PB term is actually an approximation for an external density, so

diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp
@@ -173,7 +173,7 @@ SCFEnergy FockBuilder::trace(OrbitalVector &Phi, const Nuclei &nucs) {
         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);
+        std::tuple<double, double> Er = this->Ro->getSolver()->computeEnergies(rho_el);
 
         Er_el = std::get<0>(Er);
         Er_nuc = std::get<1>(Er);

diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h
@@ -52,7 +52,7 @@ class ReactionOperator final : public RankZeroOperator {
 
     ComplexDouble trace(OrbitalVector &Phi) { return RankZeroOperator::trace(Phi); }
 
-    SCRF *getSolver() { return this->potential->getSolver(); }
+    GPESolver *getSolver() { return this->potential->getSolver(); }
     std::shared_ptr<ReactionPotential> getPotential() { return this->potential; }
     void updateMOResidual(double const err_t) { this->potential->updateMOResidual(err_t); }
 

diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h
@@ -47,7 +47,7 @@ class ReactionPotential final : public QMPotential {
     ReactionPotential(std::unique_ptr<GPESolver> gpesolver_p, std::shared_ptr<mrchem::OrbitalVector> Phi_p);
     ~ReactionPotential() override { free(NUMBER::Total); }
 
-    SCRF *getSolver() { return this->solver.get(); }
+    GPESolver *getSolver() { return this->solver.get(); }
 
     /** @brief Updates the solver.mo_residual member variable. This variable is used to set the convergence criterion in
      * the dynamic convergence method. */

diff --git a/tests/h_lpb/h.inp b/tests/h_lpb/h.inp
@@ -16,8 +16,7 @@
   "SCRF": {
     "kain": 6,
     "max_iter": 100,
-    "dynamic_thrs": false,
-    "optimizer": "potential"
+    "dynamic_thrs": false
   },
   "Cavity": {
     "spheres": "0 2.645616384 1.0 0.0 0.2"

diff --git a/tests/h_lpb/reference/h.json b/tests/h_lpb/reference/h.json
@@ -108,7 +108,6 @@
           "kain": 6,
           "kappa_out": 0.08703231499578493,
           "max_iter": 100,
-          "optimizer": "potential",
           "poisson_prec": 0.0001,
           "solver_type": "Linearized_Poisson-Boltzmann"
         },
@@ -176,7 +175,7 @@
       "charge": -1,
       "dipole_moment": {
         "dip-1": {
-          "magnitude": 8.17913080064801e-13,
+          "magnitude": 7.715572311852292e-13,
           "r_O": [
             0.0,
             0.0,
@@ -212,58 +211,58 @@
       "multiplicity": 1,
       "orbital_energies": {
         "energy": [
-          0.04928658666294078
+          -0.11654513076836014
         ],
         "occupation": [
           2.0
         ],
         "spin": [
           "p"
         ],
-        "sum_occupied": 0.09857317332588156
+        "sum_occupied": -0.23309026153672027
       },
       "scf_energy": {
-        "E_ee": 1.220280348193974,
+        "E_ee": 1.220280348193976,
         "E_eext": 0.0,
-        "E_el": -0.618896701038129,
-        "E_en": -1.9145538948639602,
-        "E_kin": 0.9258603759485293,
+        "E_el": -0.7847284165934681,
+        "E_en": -1.9145538948639653,
+        "E_kin": 0.9258603759485317,
         "E_next": 0.0,
         "E_nn": 0.0,
-        "E_nuc": 0.10754370176892118,
-        "E_tot": -0.5113529992692079,
-        "E_x": -0.15252084470528907,
-        "E_xc": -0.4885331476130991,
-        "Er_el": -0.20942953799828384,
-        "Er_nuc": 0.10754370176892118,
-        "Er_tot": -0.10188583622936287
+        "E_nuc": 0.19089592768192742,
+        "E_tot": -0.5938324889115407,
+        "E_x": -0.15252084470528926,
+        "E_xc": -0.4885331476131019,
+        "Er_el": -0.3752612535536192,
+        "Er_nuc": 0.19089592768192742,
+        "Er_tot": -0.1843653258716918
       }
     },
     "provenance": {
       "creator": "MRChem",
       "mpi_processes": 1,
-      "nthreads": 8,
+      "nthreads": 1,
       "routine": "mrchem.x",
-      "total_cores": 8,
+      "total_cores": 1,
       "version": "1.2.0-alpha"
     },
     "rsp_calculations": null,
     "scf_calculation": {
       "initial_energy": {
-        "E_ee": 1.220280348193974,
+        "E_ee": 1.220280348193976,
         "E_eext": 0.0,
-        "E_el": -0.618896701038129,
-        "E_en": -1.9145538948639602,
-        "E_kin": 0.9258603759485293,
+        "E_el": -0.7847284165934681,
+        "E_en": -1.9145538948639653,
+        "E_kin": 0.9258603759485317,
         "E_next": 0.0,
         "E_nn": 0.0,
-        "E_nuc": 0.10754370176892118,
-        "E_tot": -0.5113529992692079,
-        "E_x": -0.15252084470528907,
-        "E_xc": -0.4885331476130991,
-        "Er_el": -0.20942953799828384,
-        "Er_nuc": 0.10754370176892118,
-        "Er_tot": -0.10188583622936287
+        "E_nuc": 0.19089592768192742,
+        "E_tot": -0.5938324889115407,
+        "E_x": -0.15252084470528926,
+        "E_xc": -0.4885331476131019,
+        "Er_el": -0.3752612535536192,
+        "Er_nuc": 0.19089592768192742,
+        "Er_tot": -0.1843653258716918
       },
       "success": true
     },