Refactor the Permittivity function

src/driver.cpp

@@ -1047,7 +1047,6 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         auto poisson_prec = json_fock["reaction_operator"]["poisson_prec"];
         auto P_p = std::make_shared<PoissonOperator>(*MRA, poisson_prec);
         auto D_p = std::make_shared<mrcpp::ABGVOperator<3>>(*MRA, 0.0, 0.0);
-        auto cavity_p = mol.getCavity_p();
 
         auto kain = json_fock["reaction_operator"]["kain"];
         auto max_iter = json_fock["reaction_operator"]["max_iter"];
@@ -1057,7 +1056,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         auto eps_o = json_fock["reaction_operator"]["epsilon_out"];
         auto formulation = json_fock["reaction_operator"]["formulation"];
 
-        Permittivity dielectric_func(*cavity_p, eps_i, eps_o, formulation);
+        Permittivity dielectric_func(mol.getCavity_p(), eps_i, eps_o, formulation);
         dielectric_func.printParameters();
         Density rho_nuc(false);
         rho_nuc = chemistry::compute_nuclear_density(poisson_prec, nuclei, 100);

src/environment/CMakeLists.txt

@@ -2,4 +2,5 @@ target_sources(mrchem PRIVATE
     ${CMAKE_CURRENT_SOURCE_DIR}/Cavity.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Permittivity.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/SCRF.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/StepFunction.cpp
     )

src/environment/Permittivity.cpp

@@ -29,82 +29,20 @@
 
 namespace mrchem {
 
-Permittivity::Permittivity(const mrchem::Cavity cavity, double epsilon_in, double epsilon_out, std::string formulation)
-        : epsilon_in(epsilon_in)
-        , epsilon_out(epsilon_out)
-        , formulation(formulation)
-        , cavity(cavity) {}
+Permittivity::Permittivity(std::shared_ptr<mrchem::Cavity> cavity, double epsilon_in, double epsilon_out, std::string formulation)
+        : StepFunction(cavity, epsilon_in, epsilon_out)
+        , formulation(formulation) {}
 
 double Permittivity::evalf(const mrcpp::Coord<3> &r) const {
-    auto epsilon = epsilon_in * std::exp(std::log(epsilon_out / epsilon_in) * (1 - this->cavity.evalf(r)));
-    if (inverse) {
-        return 1 / epsilon;
-    } else {
-        return epsilon;
-    }
+    auto c_pin = this->cavity;
+    return this->in * std::exp(std::log(this->out / this->in) * (1 - c_pin->evalf(r)));
 }
 
-void Permittivity::printParameters() const {
-    // Collect relevant quantities
-    auto coords = this->cavity.getCoordinates();
-    auto radii = this->cavity.getRadii();
-    auto radii_0 = this->cavity.getOriginalRadii();
-    auto alphas = this->cavity.getRadiiScalings();
-    auto sigmas = this->cavity.getWidths();
-    auto betas = this->cavity.getWidthScalings();
-
-    // Set widths
-    auto w0 = mrcpp::Printer::getWidth() - 1;
-    auto w1 = 5;
-    auto w2 = 9;
-    auto w3 = 6;
-    auto w4 = 10;
-    auto w5 = w0 - w1 - w2 - 3 * w3 - 3 * w4;
-
-    // Build table column headers
-    std::stringstream o_head;
-    o_head << std::setw(w1) << "N";
-    o_head << std::setw(w2) << "R_0";
-    o_head << std::setw(w3 + 1) << "Alpha";
-    o_head << std::setw(w3 - 1) << "Beta";
-    o_head << std::setw(w3) << "Sigma";
-    o_head << std::setw(w5) << "Radius";
-    o_head << std::setw(w4) << "x";
-    o_head << std::setw(w4) << "y";
-    o_head << std::setw(w4) << "z";
-
-    // Print
+void Permittivity::printHeader() const {
     mrcpp::print::header(0, "Solvation Cavity");
-    print_utils::text(0, "Formulation", getFormulation(), true);
-    print_utils::scalar(0, "Dielectric constant", getEpsIn(), "(in)", 6);
-    print_utils::scalar(0, "", getEpsOut(), "(out)", 6);
-    mrcpp::print::separator(0, '-');
-    println(0, o_head.str());
-    mrcpp::print::separator(0, '-');
-    for (auto i = 0; i < coords.size(); i++) {
-        auto coord = coords[i];
-        auto x = coord[0];
-        auto y = coord[1];
-        auto z = coord[2];
-        auto r = radii[i];
-        auto r_0 = radii_0[i];
-        auto alpha = alphas[i];
-        auto beta = betas[i];
-        auto sigma = sigmas[i];
-
-        std::stringstream o_coord;
-        o_coord << std::setw(w1) << i;
-        o_coord << std::setw(w2) << std::setprecision(4) << std::fixed << r_0;
-        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << alpha;
-        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << beta;
-        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << sigma << "  ->";
-        o_coord << std::setw(w5 - 4) << std::setprecision(4) << std::fixed << r;
-        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << x;
-        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << y;
-        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << z;
-        println(0, o_coord.str());
-    }
-    mrcpp::print::separator(0, '=', 2);
+    print_utils::text(0, "Formulation", this->formulation, true);
+    print_utils::scalar(0, "Dielectric constant", getValueIn(), "(in)", 6);
+    print_utils::scalar(0, "", getValueOut(), "(out)", 6);
 }
 
 } // namespace mrchem

src/environment/Permittivity.h

@@ -26,6 +26,7 @@
 #pragma once
 
 #include "Cavity.h"
+#include "StepFunction.h"
 #include "utils/print_utils.h"
 #include <MRCPP/MWFunctions>
 #include <MRCPP/Printer>
@@ -46,7 +47,7 @@ namespace mrchem {
 
 class Cavity;
 
-class Permittivity final : public mrcpp::RepresentableFunction<3> {
+class Permittivity final : public StepFunction {
 public:
     /** @brief Standard constructor. Initializes the #cavity, #epsilon_in and #epsilon_out with the input parameters.
      *  @param cavity interlocking spheres of Cavity class.
@@ -55,7 +56,7 @@ class Permittivity final : public mrcpp::RepresentableFunction<3> {
      *  @param formulation Decides which formulation of the #Permittivity function to implement, only exponential
      * available as of now.
      */
-    Permittivity(const Cavity cavity, double epsilon_in, double epsilon_out, std::string formulation);
+    Permittivity(std::shared_ptr<Cavity> cavity, double epsilon_in, double epsilon_out, std::string formulation);
 
     /** @brief Evaluates Permittivity at a point in 3D space with respect to the state of #inverse.
      *  @param r coordinates of a 3D point in space.
@@ -64,42 +65,10 @@ class Permittivity final : public mrcpp::RepresentableFunction<3> {
      */
     double evalf(const mrcpp::Coord<3> &r) const override;
 
-    /** @brief Changes the value of #inverse. */
-    void flipFunction(bool is_inverse) { this->inverse = is_inverse; }
-
-    /** @brief Returns the current state of #inverse. */
-    auto isInverse() const { return this->inverse; }
-
-    /** @brief Calls the Cavity::getCoordinates() method of the #cavity instance. */
-    auto getCoordinates() const { return this->cavity.getCoordinates(); }
-
-    /** @brief Calls the Cavity::getRadii() method of the #cavity instance. */
-    auto getRadii() const { return this->cavity.getRadii(); }
-
-    /** @brief Calls the Cavity::getGradVector() method of the #cavity instance. */
-    auto getGradVector() const { return this->cavity.getGradVector(); }
-
-    /** @brief Returns the value of #epsilon_in. */
-    auto getEpsIn() const { return this->epsilon_in; }
-
-    /** @brief Returns the value of #epsilon_out. */
-    auto getEpsOut() const { return this->epsilon_out; }
-
-    /** @brief Returns the cavity */
-    Cavity getCavity() const { return this->cavity; }
-
-    /** @brief Returns the formulation */
-    std::string getFormulation() const { return this->formulation; }
-
-    /** @brief Print parameters */
-    void printParameters() const;
+    void printHeader() const override;
 
 private:
-    bool inverse = false;    //!< State of #evalf
-    double epsilon_in;       //!< Dielectric constant describing the permittivity of free space.
-    double epsilon_out;      //!< Dielectric constant describing the permittivity of the solvent.
-    std::string formulation; //!< Formulation of the permittivity function, only exponential is used as of now.
-    Cavity cavity;           //!< A Cavity class instance.
+    std::string formulation = "exponential";
 };
 
 } // namespace mrchem

src/environment/SCRF.cpp

@@ -96,7 +96,8 @@ void SCRF::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunctio
     resetComplexFunction(out_gamma);
 
     for (int d = 0; d < 3; d++) {
-        mrcpp::AnalyticFunction<3> d_cav(this->epsilon.getGradVector()[d]);
+        auto C_pin = this->epsilon.getCavity();
+        mrcpp::AnalyticFunction<3> d_cav(C_pin->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
@@ -107,7 +108,7 @@ void SCRF::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunctio
         }
     }
 
-    out_gamma.rescale(std::log((epsilon.getEpsIn() / epsilon.getEpsOut())) * (1.0 / (4.0 * mrcpp::pi)));
+    out_gamma.rescale(std::log((epsilon.getValueIn() / epsilon.getValueOut())) * (1.0 / (4.0 * mrcpp::pi)));
     mrcpp::clear(d_V, true);
 }
 
@@ -118,13 +119,11 @@ mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction &
     mrcpp::ComplexFunction Vr_np1;
     Vr_np1.alloc(NUMBER::Real);
 
-    this->epsilon.flipFunction(true);
-
+    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(Phi, rho_tot);
 
-    mrcpp::cplxfunc::multiply(first_term, rho_tot, this->epsilon, this->apply_prec);
-    this->epsilon.flipFunction(false);
+    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);
     rho_tot.free(NUMBER::Real);

src/environment/StepFunction.cpp

@@ -0,0 +1,107 @@
+/*
+ * MRChem, a numerical real-space code for molecular electronic structure
+ * calculations within the self-consistent field (SCF) approximations of quantum
+ * chemistry (Hartree-Fock and Density Functional Theory).
+ * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+ *
+ * This file is part of MRChem.
+ *
+ * MRChem is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MRChem is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with MRChem.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRChem, see:
+ * <https://mrchem.readthedocs.io/>
+ */
+
+#include <MRCPP/MWFunctions>
+
+#include "Cavity.h"
+#include "StepFunction.h"
+
+namespace mrchem {
+
+StepFunction::StepFunction(std::shared_ptr<mrchem::Cavity> cavity, double val_in, double val_out)
+        : in(val_in)
+        , out(val_out)
+        , cavity{std::move(cavity)} {}
+
+void StepFunction::printParameters() const {
+    // Collect relevant quantities
+    auto c_pin = this->cavity;
+
+    auto coords = c_pin->getCoordinates();
+    auto radii = c_pin->getRadii();
+    auto radii_0 = c_pin->getOriginalRadii();
+    auto alphas = c_pin->getRadiiScalings();
+    auto sigmas = c_pin->getWidths();
+    auto betas = c_pin->getWidthScalings();
+
+    // Set widths
+    auto w0 = mrcpp::Printer::getWidth() - 1;
+    auto w1 = 5;
+    auto w2 = 9;
+    auto w3 = 6;
+    auto w4 = 10;
+    auto w5 = w0 - w1 - w2 - 3 * w3 - 3 * w4;
+
+    // Build table column headers
+    std::stringstream o_head;
+    o_head << std::setw(w1) << "N";
+    o_head << std::setw(w2) << "R_0";
+    o_head << std::setw(w3 + 1) << "Alpha";
+    o_head << std::setw(w3 - 1) << "Beta";
+    o_head << std::setw(w3) << "Sigma";
+    o_head << std::setw(w5) << "Radius";
+    o_head << std::setw(w4) << "x";
+    o_head << std::setw(w4) << "y";
+    o_head << std::setw(w4) << "z";
+
+    // Print
+    printHeader();
+    mrcpp::print::separator(0, '-');
+    println(0, o_head.str());
+    mrcpp::print::separator(0, '-');
+    for (auto i = 0; i < coords.size(); i++) {
+        auto coord = coords[i];
+        auto x = coord[0];
+        auto y = coord[1];
+        auto z = coord[2];
+        auto r = radii[i];
+        auto r_0 = radii_0[i];
+        auto alpha = alphas[i];
+        auto beta = betas[i];
+        auto sigma = sigmas[i];
+
+        std::stringstream o_coord;
+        o_coord << std::setw(w1) << i;
+        o_coord << std::setw(w2) << std::setprecision(4) << std::fixed << r_0;
+        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << alpha;
+        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << beta;
+        o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << sigma << "  ->";
+        o_coord << std::setw(w5 - 4) << std::setprecision(4) << std::fixed << r;
+        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << x;
+        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << y;
+        o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << z;
+        println(0, o_coord.str());
+    }
+    mrcpp::print::separator(0, '=', 2);
+}
+
+void StepFunction::printHeader() const {
+    mrcpp::print::header(0, "Step Function of Cavity values");
+    print_utils::text(0, "Formulation", "Standard", true);
+    print_utils::scalar(0, "Value inside Cavity", getValueIn(), "(in)", 6);
+    print_utils::scalar(0, "Value outside Cavity", getValueOut(), "(out)", 6);
+}
+
+} // namespace mrchem