diff --git a/src/driver.cpp b/src/driver.cpp index cdccf7fc0..578034711 100644 --- a/src/driver.cpp +++ b/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(*MRA, poisson_prec); auto D_p = std::make_shared>(*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); diff --git a/src/environment/CMakeLists.txt b/src/environment/CMakeLists.txt index 6d855f2af..a61dc42ec 100644 --- a/src/environment/CMakeLists.txt +++ b/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 ) diff --git a/src/environment/Permittivity.cpp b/src/environment/Permittivity.cpp index fb2819765..b0cf01716 100644 --- a/src/environment/Permittivity.cpp +++ b/src/environment/Permittivity.cpp @@ -29,82 +29,13 @@ 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 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; - } -} - -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 - 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); + auto c_pin = this->cavity; + return this->in * std::exp(std::log(this->out / this->in) * (1 - c_pin->evalf(r))); } } // namespace mrchem diff --git a/src/environment/Permittivity.h b/src/environment/Permittivity.h index f7f61b44a..2d04927fd 100644 --- a/src/environment/Permittivity.h +++ b/src/environment/Permittivity.h @@ -26,6 +26,7 @@ #pragma once #include "Cavity.h" +#include "StepFunction.h" #include "utils/print_utils.h" #include #include @@ -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, 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 { detail::print_header("Solvation Cavity", this->formulation, getValueIn(), getValueOut()); } 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 diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 7ac30fc6e..8aaa1eff2 100644 --- a/src/environment/SCRF.cpp +++ b/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); diff --git a/src/environment/StepFunction.cpp b/src/environment/StepFunction.cpp new file mode 100644 index 000000000..43b628122 --- /dev/null +++ b/src/environment/StepFunction.cpp @@ -0,0 +1,108 @@ +/* + * 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ +#include "StepFunction.h" + +#include + +#include "Cavity.h" + +namespace mrchem { +namespace detail { +void print_header(const std::string &header, const std::string &formulation, double in_value, double out_value) { + mrcpp::print::header(0, header); + print_utils::text(0, "Formulation", formulation, true); + print_utils::scalar(0, "Value inside Cavity", in_value, "(in)", 6); + print_utils::scalar(0, "Value outside Cavity", out_value, "(out)", 6); +} +} // namespace detail + +StepFunction::StepFunction(std::shared_ptr 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); +} + +} // namespace mrchem diff --git a/src/environment/StepFunction.h b/src/environment/StepFunction.h new file mode 100644 index 000000000..00bfa58e9 --- /dev/null +++ b/src/environment/StepFunction.h @@ -0,0 +1,80 @@ +/* + * 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#pragma once +#include + +#include +#include + +#include "Cavity.h" +#include "utils/print_utils.h" + +namespace mrchem { + +namespace detail { +void print_header(const std::string &header, const std::string &formulation, double in_value, double out_value); +} // namespace detail + +class Cavity; + +/** @class StepFunction + * + * @brief StepFunction function related to a Cavity function. + * The StepFunction class represents the following equation + * \f[ + * S(\mathbf{r}) = \begin{cases} S_{in} & \text{if } \mathbf{r} inside C \\ + * S_{out} & \text{if } \mathbf{r} outside C + * \end{cases} + * \f] + * where \f$\mathbf{r}\f$ is the coordinate of a point in 3D space, \f$ C \f$ is the #cavity function, + * and \f$S_{in}\f$ and \f$ S_{out} \f$ are the values of the function inside and outside the #cavity respectively. + */ +class StepFunction : public mrcpp::RepresentableFunction<3> { +public: + /** @brief Standard constructor. Initializes the #cavity, #in and #out with the input parameters. + * @param cavity interlocking spheres of Cavity class. + * @param epsilon_in permittivity inside the #cavity. + * @param epsilon_out permittivity outside the #cavity. + * available as of now. + */ + StepFunction(std::shared_ptr cavity, double val_in, double val_out); + + auto getValueIn() const { return this->in; } + auto getValueOut() const { return this->out; } + + std::shared_ptr getCavity() const { return this->cavity; } + + void printParameters() const; + + virtual void printHeader() const { detail::print_header("Step function of Cavity values", "Standard", getValueIn(), getValueOut()); } + +protected: + double in; //!< Value of the function inside the #cavity. + double out; //!< value of the function outside the #cavity. + std::shared_ptr cavity; //!< A Cavity class instance. +}; + +} // namespace mrchem diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp index 045494098..d381e46c5 100644 --- a/tests/solventeffect/reaction_operator.cpp +++ b/tests/solventeffect/reaction_operator.cpp @@ -64,7 +64,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") { double slope = 0.2; std::vector radius = {1.0}; std::vector> coords = {{0.0, 0.0, 0.0}}; - Cavity sphere(coords, radius, slope); + auto sphere = std::make_shared(coords, radius, slope); // initialize dielectric function double eps_in = 1.0;