From 2eab5490cce7ae539c0351049bf4e4bdc25b5fd8 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Wed, 21 Mar 2018 20:16:34 -0400 Subject: [PATCH 1/3] Simplify RadialSolution --- CHANGELOG.md | 3 + src/green/InterfacesImpl.hpp | 325 +++++++++++---------------------- src/green/SphericalDiffuse.cpp | 14 +- src/green/SphericalDiffuse.hpp | 10 +- 4 files changed, 117 insertions(+), 235 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 941e187ee..296dd3a73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,6 +31,9 @@ - The `ENABLE_Fortran_API` configuration option has been renamed `TEST_Fortran_API`, since it now only triggers compilation of the `Fortran_host` test case. +- The implementation of the `RadialSolution` for the second order ODE + associated with the spherical diffuse Green's function is less heavily + `template`-d. ## [Version 1.2.0-rc1] - 2018-03-02 diff --git a/src/green/InterfacesImpl.hpp b/src/green/InterfacesImpl.hpp index 75bc56e51..aa92e65cb 100644 --- a/src/green/InterfacesImpl.hpp +++ b/src/green/InterfacesImpl.hpp @@ -38,16 +38,14 @@ #include #include "utils/MathUtils.hpp" +#include "utils/RungeKutta4.hpp" + +/*! \file InterfacesImpl.hpp */ namespace pcm { namespace green { namespace detail { -/*! \typedef StateType - * \brief state vector for the differential equation integrator - */ -typedef std::vector StateType; - /*! \typedef ProfileEvaluator * \brief sort of a function pointer to the dielectric profile evaluation function */ @@ -75,16 +73,18 @@ struct IntegratorParameters { * Provides a handle to the system of differential equations for the integrator. * The dielectric profile comes in as a boost::function object. */ -class LnTransformedRadial { +class LnTransformedRadial __final : public pcm::utils::detail::ODESystem<2> { +public: + /*! Type of the state vector of the ODE */ + typedef pcm::utils::detail::ODESystem<2>::StateType StateType; + /*! Constructor from profile evaluator and angular momentum */ + LnTransformedRadial(const ProfileEvaluator & e, int lval) : eval_(e), l_(lval) {} + private: /*! Dielectric profile function and derivative evaluation */ ProfileEvaluator eval_; /*! Angular momentum */ int l_; - -public: - /*! Constructor from profile evaluator and angular momentum */ - LnTransformedRadial(const ProfileEvaluator & e, int lval) : eval_(e), l_(lval) {} /*! Provides a functor for the evaluation of the system * of first-order ODEs needed by Boost.Odeint * The second-order ODE and the system of first-order ODEs @@ -93,7 +93,7 @@ class LnTransformedRadial { * \param[out] drhodr state vector holding the first and second derivative * \param[in] y logarithmic position on the integration grid */ - void operator()(const StateType & rho, StateType & drhodr, const double y) { + virtual void RHS(const StateType & rho, StateType & drhodr, const double y) const { // Evaluate the dielectric profile double eps = 0.0, epsPrime = 0.0; pcm::tie(eps, epsPrime) = eval_(std::exp(y)); @@ -107,76 +107,49 @@ class LnTransformedRadial { }; } // namespace detail -using detail::ProfileEvaluator; using detail::IntegratorParameters; +using detail::ProfileEvaluator; -/*! \file InterfacesImpl.hpp - * \class RadialFunction +/*! \class RadialFunction * \brief represents solutions to the radial 2nd order ODE * \author Roberto Di Remigio - * \date 2015 - * \tparam StateVariable type of the state variable used in the ODE solver - * \tparam ODESystem system of 1st order ODEs replacing the 2nd order ODE - * \tparam IndependentSolution encodes which type of radial solution + * \date 2018 + * \tparam ODE system of 1st order ODEs replacing the 2nd order ODE + * + * The sign of the integration interval, i.e. the sign of the difference + * between the minimum and the maximum of the interval, is used to discriminate + * between the two independent solutions (zeta-type and omega-type) of the + * differential equation: + * - When the sign is positive (y_min_ < y_max_), we are integrating **forwards** + * for the zeta-type solution. + * - When the sign is negative (y_max_ < y_min_), we are integrating **backwards** + * for the omega-type solution. + * + * This is based on an idea Luca had in Wuhan/Chongqing for the planar + * diffuse interface. + * With respect to the previous, much more heavily template-d implementation, + * it introduces a bit more if-else if-else branching in the function_impl and + * derivative_impl functions. + * This is largely offset by the better readability and reduced code duplication. */ -template class IndependentSolution> -class RadialFunction __final { +template class RadialFunction __final { public: - RadialFunction() : solution_(IndependentSolution()) {} + RadialFunction() : L_(0), y_min_(0.0), y_max_(0.0), y_sign_(1) {} RadialFunction(int l, - double y0, - double yinf, + double ymin, + double ymax, const ProfileEvaluator & eval, const IntegratorParameters & parms) - : solution_(IndependentSolution(l, - y0, - yinf, - eval, - parms)) {} - ~RadialFunction() {} - /*! \brief Returns value of function and its first derivative at given point - * \param[in] point evaluation point - */ - pcm::tuple operator()(double point) const { - return solution_(point); - } - friend std::ostream & operator<<(std::ostream & os, RadialFunction & obj) { - os << obj.solution_; - return os; - } - -private: - /// Independent solution to the radial equation - IndependentSolution solution_; -}; - -/*! \file InterfacesImpl.hpp - * \class Zeta - * \brief 1st solution to the radial second order ODE, with r^l behaviour in the - * origin - * \author Roberto Di Remigio - * \date 2015 - * \tparam StateVariable type of the state variable used in the ODE solver - * \tparam ODESystem system of 1st order ODEs replacing the 2nd order ODE - */ -template class Zeta __final { -public: - Zeta() : L_(0), y_0_(0.0), y_infinity_(0.0) {} - Zeta(int l, - double y0, - double yinf, - const ProfileEvaluator & eval, - const IntegratorParameters & parms) - : L_(l), y_0_(y0), y_infinity_(yinf) { + : L_(l), + y_min_(ymin), + y_max_(ymax), + y_sign_(pcm::utils::sign(y_max_ - y_min_)) { compute(eval, parms); } - ~Zeta() {} pcm::tuple operator()(double point) const { return pcm::make_tuple(function_impl(point), derivative_impl(point)); } - friend std::ostream & operator<<(std::ostream & os, Zeta & obj) { + friend std::ostream & operator<<(std::ostream & os, RadialFunction & obj) { for (size_t i = 0; i < obj.function_[0].size(); ++i) { os << obj.function_[0][i] << " " << obj.function_[1][i] << " " << obj.function_[2][i] << std::endl; @@ -185,205 +158,113 @@ template class Zeta __final { } private: - typedef pcm::array RadialSolution; - /// Angular momentum of the solution + typedef typename ODE::StateType StateType; + /*! Angular momentum of the solution */ int L_; - /// Lower bound of the integration interval - double y_0_; - /// Upper bound of the integration interval - double y_infinity_; - /// The actual data: grid, function value and first derivative values - RadialSolution function_; + /*! Lower bound of the integration interval */ + double y_min_; + /*! Upper bound of the integration interval */ + double y_max_; + /*! The sign of the integration interval */ + double y_sign_; + /*! The actual data: grid, function value and first derivative values */ + pcm::array, 3> function_; /*! Reports progress of differential equation integrator */ - void push_back(const StateVariable & x, double y) { + void push_back(const StateType & x, double y) { function_[0].push_back(y); function_[1].push_back(x[0]); function_[2].push_back(x[1]); } - /*! \brief Calculates 1st radial solution, i.e. the one with r^l behavior + /*! \brief Calculates radial solution * \param[in] eval dielectric profile evaluator function object * \param[in] parms parameters for the integrator - */ - void compute(const ProfileEvaluator & eval, const IntegratorParameters & parms) { - namespace odeint = boost::numeric::odeint; - odeint::runge_kutta4 stepper; - - ODESystem system(eval, L_); - // Holds the initial conditions - StateVariable init_zeta(2); - // Set initial conditions - init_zeta[0] = L_ * y_0_; - init_zeta[1] = L_; - odeint::integrate_const( - stepper, - system, - init_zeta, - y_0_, - y_infinity_, - parms.observer_step_, - pcm::bind( - &Zeta::push_back, this, pcm::_1, pcm::_2)); - } - /*! \brief Returns value of function at given point - * \param[in] point evaluation point * - * We first check if point is below y_0_, if yes we use - * the asymptotic form L*y in point. - */ - double function_impl(double point) const { - double zeta = 0.0; - if (point <= y_0_) { - zeta = L_ * point; - } else { - zeta = utils::splineInterpolation(point, function_[0], function_[1]); - } - return zeta; - } - /*! \brief Returns value of 1st derivative of function at given point - * \param[in] point evaluation point - * - * Below y_0_, the as asymptotic form L is used. Otherwise we interpolate. - */ - double derivative_impl(double point) const { - double zeta = 0.0; - if (point <= y_0_) { - zeta = L_; - } else { - zeta = utils::splineInterpolation(point, function_[0], function_[2]); - } - return zeta; - } -}; - -/*! \file InterfacesImpl.hpp - * \class Omega - * \brief 2nd solution to the radial second order ODE, with r^(-l-1) behaviour at - * infinity - * \author Roberto Di Remigio - * \date 2015 - * \tparam StateVariable type of the state variable used in the ODE solver - * \tparam ODESystem system of 1st order ODEs replacing the 2nd order ODE - */ -template class Omega __final { -public: - Omega() : L_(0), y_0_(0.0), y_infinity_(0.0) {} - Omega(int l, - double y0, - double yinf, - const ProfileEvaluator & eval, - const IntegratorParameters & parms) - : L_(l), y_0_(y0), y_infinity_(yinf) { - compute(eval, parms); - } - ~Omega() {} - pcm::tuple operator()(double point) const { - return pcm::make_tuple(function_impl(point), derivative_impl(point)); - } - friend std::ostream & operator<<(std::ostream & os, Omega & obj) { - for (size_t i = 0; i < obj.function_[0].size(); ++i) { - os << obj.function_[0][i] << " " << obj.function_[1][i] << " " - << obj.function_[2][i] << std::endl; - } - return os; - } - -private: - typedef pcm::array RadialSolution; - /// Angular momentum of the solution - int L_; - /// Lower bound of the integration interval - double y_0_; - /// Upper bound of the integration interval - double y_infinity_; - /// The actual data: grid, function value and first derivative values - RadialSolution function_; - /*! Reports progress of differential equation integrator */ - void push_back(const StateVariable & x, double y) { - function_[0].push_back(y); - function_[1].push_back(x[0]); - function_[2].push_back(x[1]); - } - /*! \brief calculates 2nd radial solution, i.e. the one with r^(-l-1) behavior - * \param[in] eval dielectric profile evaluator function object - * \param[in] parms parameters for the integrator + * This function discriminates between the first, i.e. the one with r^l + * behavior, and the second radial solution, i.e. the one with r^(-l-1) + * behavior, based on the sign of the integration interval y_sign_. */ void compute(const ProfileEvaluator & eval, const IntegratorParameters & parms) { namespace odeint = boost::numeric::odeint; - odeint::runge_kutta4 stepper; + odeint::runge_kutta4 stepper; - ODESystem system(eval, L_); + ODE system(eval, L_); // Holds the initial conditions - StateVariable init_omega(2); + StateType init; // Set initial conditions - init_omega[0] = -(L_ + 1) * y_infinity_; - init_omega[1] = -(L_ + 1); - // Notice that we integrate BACKWARDS, so we pass -step to integrate_adaptive + if (y_sign_ > 0.0) { // zeta-type solution + init[0] = y_sign_ * L_ * y_min_; + init[1] = y_sign_ * L_; + } else { // omega-type solution + init[0] = y_sign_ * (L_ + 1) * y_min_; + init[1] = y_sign_ * (L_ + 1); + } odeint::integrate_const( stepper, system, - init_omega, - y_infinity_, - y_0_, - -parms.observer_step_, - pcm::bind( - &Omega::push_back, this, pcm::_1, pcm::_2)); -// Reverse order of StateVariable-s in RadialSolution -// this ensures that they are in ascending order, as later expected by -// function_impl and derivative_impl + init, + y_min_, + y_max_, + y_sign_ * parms.observer_step_, + pcm::bind(&RadialFunction::push_back, this, pcm::_1, pcm::_2)); + // clang-format off + // Reverse order of function_ if omega-type solution was computed + // this ensures that they are in ascending order, as later expected by + // function_impl and derivative_impl + if (y_sign_ < 0.0) { #ifdef HAS_CXX11 - for (auto & comp : function_) { + for (auto & comp : function_) { #else /* HAS_CXX11 */ - BOOST_FOREACH (StateVariable & comp, function_) { + BOOST_FOREACH (std::vector & comp, function_) { #endif /* HAS_CXX11 */ - std::reverse(comp.begin(), comp.end()); + std::reverse(comp.begin(), comp.end()); + } + } } - } + // clang-format on /*! \brief Returns value of function at given point * \param[in] point evaluation point * - * We first check if point is above r_infinity_, if yes we use - * the asymptotic form -(L+1)*y in point. + * We first check if point is below y_min_, if yes we use + * the asymptotic form. */ double function_impl(double point) const { - double omega = 0.0; - if (point >= y_infinity_) { - omega = -(L_ + 1) * point; + double val = 0.0; + if (point <= y_min_ && y_sign_ > 0.0) { // Asymptotic zeta-type solution + val = y_sign_ * L_ * point; + } else if (point >= y_min_ && y_sign_ < 0.0) { // Asymptotic omega-type solution + val = y_sign_ * (L_ + 1) * point; } else { - omega = utils::splineInterpolation(point, function_[0], function_[1]); + val = utils::splineInterpolation(point, function_[0], function_[1]); } - return omega; + return val; } /*! \brief Returns value of 1st derivative of function at given point * \param[in] point evaluation point * - * We first check if point is above r_infinity_, if yes we use - * the asymptotic form -(L+1) in point. + * Below y_min_, the asymptotic form is used. Otherwise we interpolate. */ double derivative_impl(double point) const { - double omega = 0.0; - if (point >= y_infinity_) { - omega = -(L_ + 1); + double val = 0.0; + if (point <= y_min_ && y_sign_ > 0.0) { // Asymptotic zeta-type solution + val = y_sign_ * L_; + } else if (point >= y_min_ && y_sign_ < 0.0) { // Asymptotic omega-type solution + val = y_sign_ * (L_ + 1); } else { - omega = utils::splineInterpolation(point, function_[0], function_[2]); + val = utils::splineInterpolation(point, function_[0], function_[2]); } - return omega; + return val; } }; /*! \brief Write contents of a RadialFunction to file - * \param[in] f RadialSolution whose contents have to be printed + * \param[in] f solution to the radial 2nd order ODE to print * \param[in] fname name of the file * \author Roberto Di Remigio - * \date 2015 - * \tparam StateVariable type of the state variable used in the ODE solver - * \tparam ODESystem system of 1st order ODEs replacing the 2nd order ODE - * \tparam IndependentSolution encodes which type of radial solution + * \date 2018 + * \tparam ODE system of 1st order ODEs replacing the 2nd order ODE */ -template class IndependentSolution> -void writeToFile(RadialFunction & f, - const std::string & fname) { +template +void writeToFile(RadialFunction & f, const std::string & fname) { std::ofstream fout; fout.open(fname.c_str()); fout << f << std::endl; diff --git a/src/green/SphericalDiffuse.cpp b/src/green/SphericalDiffuse.cpp index 9480eded2..ab0de288e 100644 --- a/src/green/SphericalDiffuse.cpp +++ b/src/green/SphericalDiffuse.cpp @@ -254,20 +254,20 @@ void SphericalDiffuse::initSphericalDiffuse() { ProfileEvaluator eval_ = pcm::bind(&ProfilePolicy::operator(), this->profile_, pcm::_1); - zetaC_ = RadialFunction( - maxLC_, y_0_, y_infinity_, eval_, params_); - omegaC_ = RadialFunction( - maxLC_, y_0_, y_infinity_, eval_, params_); + zetaC_ = + RadialFunction(maxLC_, y_0_, y_infinity_, eval_, params_); + omegaC_ = + RadialFunction(maxLC_, y_infinity_, y_0_, eval_, params_); zeta_.reserve(maxLGreen_ + 1); omega_.reserve(maxLGreen_ + 1); for (int L = 0; L <= maxLGreen_; ++L) { - RadialFunction tmp_zeta_( + RadialFunction tmp_zeta_( L, y_0_, y_infinity_, eval_, params_); zeta_.push_back(tmp_zeta_); } for (int L = 0; L <= maxLGreen_; ++L) { - RadialFunction tmp_omega_( - L, y_0_, y_infinity_, eval_, params_); + RadialFunction tmp_omega_( + L, y_infinity_, y_0_, eval_, params_); omega_.push_back(tmp_omega_); } } diff --git a/src/green/SphericalDiffuse.hpp b/src/green/SphericalDiffuse.hpp index b90ea826e..de7914190 100644 --- a/src/green/SphericalDiffuse.hpp +++ b/src/green/SphericalDiffuse.hpp @@ -213,14 +213,12 @@ class SphericalDiffuse __final : public GreensFunction { /*! \brief First independent radial solution, used to build Green's function. * \note The vector has dimension maxLGreen_ and has r^l behavior */ - std::vector > - zeta_; + std::vector > zeta_; /*! \brief Second independent radial solution, used to build Green's function. * \note The vector has dimension maxLGreen_ and has r^(-l-1) behavior */ - std::vector > - omega_; + std::vector > omega_; /*! \brief Returns L-th component of the radial part of the Green's function * \param[in] L angular momentum @@ -246,13 +244,13 @@ class SphericalDiffuse __final : public GreensFunction { /*! \brief First independent radial solution, used to build coefficient. * \note This is needed to separate the Coulomb singularity and has r^l behavior */ - RadialFunction zetaC_; + RadialFunction zetaC_; /*! \brief Second independent radial solution, used to build coefficient. * \note This is needed to separate the Coulomb singularity and has r^(-l-1) * behavior */ - RadialFunction omegaC_; + RadialFunction omegaC_; /*! \brief Returns coefficient for the separation of the Coulomb singularity * \param[in] sp first point From 174b0271b320aec3617f9033f9ae6acdd9148b52 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Thu, 22 Mar 2018 00:54:50 -0400 Subject: [PATCH 2/3] Replace Odeint with in-house RK4 --- .gitattributes | 2 +- CHANGELOG.md | 1 + src/green/InterfacesImpl.hpp | 97 +++++++------- src/green/SphericalDiffuse.cpp | 35 +++-- src/interface/Meddle.cpp | 2 +- src/utils/CMakeLists.txt | 1 + src/utils/MathUtils.hpp | 227 +-------------------------------- src/utils/Molecule.cpp | 6 +- src/utils/Molecule.hpp | 3 +- src/utils/RungeKutta.hpp | 197 ++++++++++++++++++++++++++++ src/utils/RungeKutta4.cpp | 134 ------------------- src/utils/RungeKutta4.hpp | 82 ------------ src/utils/SplineFunction.hpp | 49 ++++++- src/utils/Symmetry.cpp | 100 ++++++++++++--- src/utils/Symmetry.hpp | 166 ++++++++++++++++++------ tests/TestingMolecules.hpp | 54 ++++---- tests/green/CMakeLists.txt | 6 +- 17 files changed, 568 insertions(+), 594 deletions(-) create mode 100644 src/utils/RungeKutta.hpp delete mode 100644 src/utils/RungeKutta4.cpp delete mode 100644 src/utils/RungeKutta4.hpp diff --git a/.gitattributes b/.gitattributes index 24ab84ebc..4291ccfbb 100644 --- a/.gitattributes +++ b/.gitattributes @@ -51,7 +51,7 @@ src/utils/*.cpp licensefile=.githooks/LICENSE-C++ src/utils/cnpy.hpp !licensefile src/utils/cnpy.cpp !licensefile src/utils/Legendre.hpp !licensefile -src/utils/RungeKutta4.* !licensefile +src/utils/RungeKutta.hpp !licensefile # tests tests/unit_tests.cpp licensefile=.githooks/LICENSE-C++ diff --git a/CHANGELOG.md b/CHANGELOG.md index 296dd3a73..678662f10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,6 +34,7 @@ - The implementation of the `RadialSolution` for the second order ODE associated with the spherical diffuse Green's function is less heavily `template`-d. +- The dependency on Boost.Odeint has been dropped. ## [Version 1.2.0-rc1] - 2018-03-02 diff --git a/src/green/InterfacesImpl.hpp b/src/green/InterfacesImpl.hpp index aa92e65cb..0a5905e10 100644 --- a/src/green/InterfacesImpl.hpp +++ b/src/green/InterfacesImpl.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include "Config.hpp" @@ -34,37 +35,39 @@ #ifndef HAS_CXX11 #include #endif -// Boost.Odeint includes -#include #include "utils/MathUtils.hpp" -#include "utils/RungeKutta4.hpp" +#include "utils/RungeKutta.hpp" +#include "utils/SplineFunction.hpp" /*! \file InterfacesImpl.hpp */ namespace pcm { namespace green { namespace detail { +/*! \brief Abstract class for an system of ordinary differential equations + * \tparam Order The order of the ordinary differential equation + * \author Roberto Di Remigio + * \date 2018 + */ +template class ODESystem { +public: + typedef pcm::array StateType; + size_t ODEorder() const { return Order; } + void operator()(const StateType & f, StateType & dfdx, const double t) const { + RHS(f, dfdx, t); + } + virtual ~ODESystem() {} + +private: + virtual void RHS(const StateType & f, StateType & dfdx, const double t) const = 0; +}; /*! \typedef ProfileEvaluator * \brief sort of a function pointer to the dielectric profile evaluation function */ typedef pcm::function(const double)> ProfileEvaluator; -/*! \struct IntegratorParameters - * \brief holds parameters for the integrator - */ -struct IntegratorParameters { - /*! Lower bound of the integration interval */ - double r_0_; - /*! Upper bound of the integration interval */ - double r_infinity_; - /*! Time step between observer calls */ - double observer_step_; - IntegratorParameters(double r0, double rinf, double step) - : r_0_(r0), r_infinity_(rinf), observer_step_(step) {} -}; - /*! \class LnTransformedRadial * \brief system of ln-transformed first-order radial differential equations * \author Roberto Di Remigio @@ -73,25 +76,25 @@ struct IntegratorParameters { * Provides a handle to the system of differential equations for the integrator. * The dielectric profile comes in as a boost::function object. */ -class LnTransformedRadial __final : public pcm::utils::detail::ODESystem<2> { +class LnTransformedRadial __final : public pcm::green::detail::ODESystem<2> { public: /*! Type of the state vector of the ODE */ - typedef pcm::utils::detail::ODESystem<2>::StateType StateType; + typedef pcm::green::detail::ODESystem<2>::StateType StateType; /*! Constructor from profile evaluator and angular momentum */ LnTransformedRadial(const ProfileEvaluator & e, int lval) : eval_(e), l_(lval) {} private: /*! Dielectric profile function and derivative evaluation */ - ProfileEvaluator eval_; + const ProfileEvaluator eval_; /*! Angular momentum */ - int l_; - /*! Provides a functor for the evaluation of the system - * of first-order ODEs needed by Boost.Odeint - * The second-order ODE and the system of first-order ODEs - * are reported in the manuscript. + const int l_; + /*! \brief Provides a functor for the evaluation of the system of first-order ODEs. * \param[in] rho state vector holding the function and its first derivative * \param[out] drhodr state vector holding the first and second derivative * \param[in] y logarithmic position on the integration grid + * + * The second-order ODE and the system of first-order ODEs + * are reported in the manuscript. */ virtual void RHS(const StateType & rho, StateType & drhodr, const double y) const { // Evaluate the dielectric profile @@ -107,7 +110,6 @@ class LnTransformedRadial __final : public pcm::utils::detail::ODESystem<2> { }; } // namespace detail -using detail::IntegratorParameters; using detail::ProfileEvaluator; /*! \class RadialFunction @@ -138,21 +140,25 @@ template class RadialFunction __final { RadialFunction(int l, double ymin, double ymax, - const ProfileEvaluator & eval, - const IntegratorParameters & parms) + double ystep, + const ProfileEvaluator & eval) : L_(l), y_min_(ymin), y_max_(ymax), y_sign_(pcm::utils::sign(y_max_ - y_min_)) { - compute(eval, parms); + compute(ystep, eval); } pcm::tuple operator()(double point) const { return pcm::make_tuple(function_impl(point), derivative_impl(point)); } friend std::ostream & operator<<(std::ostream & os, RadialFunction & obj) { for (size_t i = 0; i < obj.function_[0].size(); ++i) { - os << obj.function_[0][i] << " " << obj.function_[1][i] << " " + // clang-format off + os << std::fixed << std::left << std::setprecision(14) + << obj.function_[0][i] << " " + << obj.function_[1][i] << " " << obj.function_[2][i] << std::endl; + // clang-format on } return os; } @@ -176,21 +182,19 @@ template class RadialFunction __final { function_[2].push_back(x[1]); } /*! \brief Calculates radial solution - * \param[in] eval dielectric profile evaluator function object - * \param[in] parms parameters for the integrator + * \param[in] step ODE integrator step + * \param[in] eval dielectric profile evaluator function object + * \return the number of integration steps * - * This function discriminates between the first, i.e. the one with r^l - * behavior, and the second radial solution, i.e. the one with r^(-l-1) - * behavior, based on the sign of the integration interval y_sign_. + * This function discriminates between the first (zeta-type), i.e. the one + * with r^l behavior, and the second (omega-type) radial solution, i.e. the + * one with r^(-l-1) behavior, based on the sign of the integration interval + * y_sign_. */ - void compute(const ProfileEvaluator & eval, const IntegratorParameters & parms) { - namespace odeint = boost::numeric::odeint; - odeint::runge_kutta4 stepper; - + size_t compute(const double step, const ProfileEvaluator & eval) { ODE system(eval, L_); - // Holds the initial conditions - StateType init; // Set initial conditions + StateType init; if (y_sign_ > 0.0) { // zeta-type solution init[0] = y_sign_ * L_ * y_min_; init[1] = y_sign_ * L_; @@ -198,14 +202,16 @@ template class RadialFunction __final { init[0] = y_sign_ * (L_ + 1) * y_min_; init[1] = y_sign_ * (L_ + 1); } - odeint::integrate_const( + pcm::utils::RungeKutta4 stepper; + size_t nSteps = pcm::utils::integrate_const( stepper, system, init, y_min_, y_max_, - y_sign_ * parms.observer_step_, + y_sign_ * step, pcm::bind(&RadialFunction::push_back, this, pcm::_1, pcm::_2)); + // clang-format off // Reverse order of function_ if omega-type solution was computed // this ensures that they are in ascending order, as later expected by @@ -218,9 +224,10 @@ template class RadialFunction __final { #endif /* HAS_CXX11 */ std::reverse(comp.begin(), comp.end()); } - } + // clang-format on } - // clang-format on + return nSteps; + } /*! \brief Returns value of function at given point * \param[in] point evaluation point * diff --git a/src/green/SphericalDiffuse.cpp b/src/green/SphericalDiffuse.cpp index ab0de288e..968af29d7 100644 --- a/src/green/SphericalDiffuse.cpp +++ b/src/green/SphericalDiffuse.cpp @@ -129,8 +129,10 @@ void SphericalDiffuse::toFile(const std::string & prefix) { writeToFile(zetaC_, tmp + "zetaC.dat"); writeToFile(omegaC_, tmp + "omegaC.dat"); for (int L = 1; L <= maxLGreen_; ++L) { - writeToFile(zeta_[L], tmp + "zeta_" + pcm::to_string(L) + ".dat"); - writeToFile(omega_[L], tmp + "omega_" + pcm::to_string(L) + ".dat"); + std::ostringstream Llabel; + Llabel << std::setw(4) << std::setfill('0') << L; + writeToFile(zeta_[L], tmp + "zeta_" + Llabel.str() + ".dat"); + writeToFile(omega_[L], tmp + "omega_" + Llabel.str() + ".dat"); } } @@ -240,35 +242,30 @@ void SphericalDiffuse::initSphericalDiffuse() { using namespace detail; // Parameters for the numerical solution of the radial differential equation - double r_0_ = 0.1; /*! Lower bound of the integration interval */ - double r_infinity_ = this->profile_.upperLimit() + - 30.0; /*! Upper bound of the integration interval */ - + /*! Lower bound of the integration interval */ + double r_0_ = 0.1; double y_0_ = std::log(r_0_); + /*! Upper bound of the integration interval */ + double r_infinity_ = this->profile_.upperLimit() + 30.0; double y_infinity_ = std::log(r_infinity_); double relative_width = this->profile_.relativeWidth(); - double observer_step_ = - 1.0e-2 * relative_width; /*! Time step between observer calls */ + /*! Time step between observer calls */ + double step_ = 1.0e-2 * relative_width; - IntegratorParameters params_(y_0_, y_infinity_, observer_step_); ProfileEvaluator eval_ = pcm::bind(&ProfilePolicy::operator(), this->profile_, pcm::_1); zetaC_ = - RadialFunction(maxLC_, y_0_, y_infinity_, eval_, params_); + RadialFunction(maxLC_, y_0_, y_infinity_, step_, eval_); omegaC_ = - RadialFunction(maxLC_, y_infinity_, y_0_, eval_, params_); + RadialFunction(maxLC_, y_infinity_, y_0_, step_, eval_); zeta_.reserve(maxLGreen_ + 1); omega_.reserve(maxLGreen_ + 1); for (int L = 0; L <= maxLGreen_; ++L) { - RadialFunction tmp_zeta_( - L, y_0_, y_infinity_, eval_, params_); - zeta_.push_back(tmp_zeta_); - } - for (int L = 0; L <= maxLGreen_; ++L) { - RadialFunction tmp_omega_( - L, y_infinity_, y_0_, eval_, params_); - omega_.push_back(tmp_omega_); + zeta_.push_back( + RadialFunction(L, y_0_, y_infinity_, step_, eval_)); + omega_.push_back( + RadialFunction(L, y_infinity_, y_0_, step_, eval_)); } } diff --git a/src/interface/Meddle.cpp b/src/interface/Meddle.cpp index 2f56d6a74..c3b03582f 100644 --- a/src/interface/Meddle.cpp +++ b/src/interface/Meddle.cpp @@ -459,7 +459,7 @@ void Meddle::initInput(int nr_nuclei, Eigen::Matrix3Xd centers = Eigen::Map(coordinates, 3, nr_nuclei); if (input_.mode() != "EXPLICIT") { - Symmetry pg = buildGroup( + Symmetry pg( symmetry_info[0], symmetry_info[1], symmetry_info[2], symmetry_info[3]); input_.molecule(detail::initMolecule(input_, pg, nr_nuclei, chg, centers)); } diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 7f9ed9784..0b530f6c0 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -9,6 +9,7 @@ LoggerImpl.hpp MathUtils.hpp Molecule.hpp QuadratureRules.hpp +RungeKutta.hpp Solvent.hpp Sphere.hpp SplineFunction.hpp diff --git a/src/utils/MathUtils.hpp b/src/utils/MathUtils.hpp index 770def5eb..9216120b9 100644 --- a/src/utils/MathUtils.hpp +++ b/src/utils/MathUtils.hpp @@ -23,70 +23,21 @@ #pragma once -#include -#include #include #include #include -#include -#include -#include #include "Config.hpp" #include #include -#include "SplineFunction.hpp" -#include "Symmetry.hpp" #include "cnpy.hpp" /*! \file MathUtils.hpp */ namespace pcm { namespace utils { -/*! \fn inline int parity(std::bitset bitrep) - * \param[in] bitrep a bitset - * \tparam nBits lenght of the input bitset - * - * Calculate the parity of the bitset as defined by: - * bitrep[0] XOR bitrep[1] XOR ... XOR bitrep[nBits-1] - */ -template inline int parity(std::bitset bitrep) { - int parity = 0; - for (size_t i = 0; i < bitrep.size(); ++i) { - parity ^= bitrep[i]; - } - return parity; -} - -/*! \fn inline double parity(unsigned int i) - * \param[in] i an integer, usually an index for an irrep or a symmetry operation - * - * Returns parity of input integer. - * The parity is defined as the result of using XOR on the bitrep - * of the given integer. For example: - * 2 -> 010 -> 0^1^0 = 1 -> -1.0 - * 6 -> 110 -> 1^1^0 = 0 -> 1.0 - * - * It can also be interpreted as the action of a given operation on the - * Cartesian axes: - * zyx Parity - * 0 000 E 1.0 - * 1 001 Oyz -1.0 - * 2 010 Oxz -1.0 - * 3 011 C2z 1.0 - * 4 100 Oxy -1.0 - * 5 101 C2y 1.0 - * 6 110 C2x 1.0 - * 7 111 i -1.0 - */ -inline double parity(unsigned int i) { - // Use a ternary if construct. If the bitset is odd return -1.0 Return +1.0 - // otherwise. - return (parity(std::bitset<3>(i)) ? -1.0 : 1.0); -} - /*! \fn inline bool isZero(double value, double threshold) * \param[in] value the value to be checked * \param[in] threshold the threshold @@ -113,110 +64,7 @@ inline bool numericalZero(double value) { return (isZero(value, 1.0e-14)); } */ template inline int sign(T val) { return (T(0) < val) - (val < T(0)); } -/*! \fn inline void symmetryBlocking(Eigen::MatrixXd & matrix, int cavitySize, int - * ntsirr, int nr_irrep) - * \param[out] matrix the matrix to be block-diagonalized - * \param[in] cavitySize the size of the cavity (size of the matrix) - * \param[in] ntsirr the size of the irreducible portion of the cavity (size of - * the blocks) - * \param[in] nr_irrep the number of irreducible representations (number of - * blocks) - */ -inline void symmetryBlocking(Eigen::MatrixXd & matrix, - PCMSolverIndex cavitySize, - PCMSolverIndex ntsirr, - int nr_irrep) { - // This function implements the simmetry-blocking of the PCM - // matrix due to point group symmetry as reported in: - // L. Frediani, R. Cammi, C. S. Pomelli, J. Tomasi and K. Ruud, J. Comput.Chem. 25, - // 375 (2003) - // u is the character table for the group (t in the paper) - Eigen::MatrixXd u = Eigen::MatrixXd::Zero(nr_irrep, nr_irrep); - for (int i = 0; i < nr_irrep; ++i) { - for (int j = 0; j < nr_irrep; ++j) { - u(i, j) = parity(i & j); - } - } - // Naming of indices: - // a, b, c, d run over the total size of the cavity (N) - // i, j, k, l run over the number of irreps (n) - // p, q, r, s run over the irreducible size of the cavity (N/n) - // Instead of forming U (T in the paper) and then perform the dense - // multiplication, we multiply block-by-block using just the u matrix. - // matrix = U * matrix * Ut; U * Ut = Ut * U = id - // First half-transformation, i.e. first_half = matrix * Ut - Eigen::MatrixXd first_half = Eigen::MatrixXd::Zero(cavitySize, cavitySize); - for (int i = 0; i < nr_irrep; ++i) { - int ioff = i * ntsirr; - for (int k = 0; k < nr_irrep; ++k) { - int koff = k * ntsirr; - for (int j = 0; j < nr_irrep; ++j) { - int joff = j * ntsirr; - double ujk = u(j, k) / nr_irrep; - for (int p = 0; p < ntsirr; ++p) { - int a = ioff + p; - for (int q = 0; q < ntsirr; ++q) { - int b = joff + q; - int c = koff + q; - first_half(a, c) += matrix(a, b) * ujk; - } - } - } - } - } - // Second half-transformation, i.e. matrix = U * first_half - matrix.setZero(cavitySize, cavitySize); - for (int i = 0; i < nr_irrep; ++i) { - int ioff = i * ntsirr; - for (int k = 0; k < nr_irrep; ++k) { - int koff = k * ntsirr; - for (int j = 0; j < nr_irrep; ++j) { - int joff = j * ntsirr; - double uij = u(i, j); - for (int p = 0; p < ntsirr; ++p) { - int a = ioff + p; - int b = joff + p; - for (int q = 0; q < ntsirr; ++q) { - int c = koff + q; - matrix(a, c) += uij * first_half(b, c); - } - } - } - } - } - // Traverse the matrix and discard numerical zeros - for (PCMSolverIndex a = 0; a < cavitySize; ++a) { - for (PCMSolverIndex b = 0; b < cavitySize; ++b) { - if (numericalZero(matrix(a, b))) { - matrix(a, b) = 0.0; - } - } - } -} - -/*! \fn inline void symmetryPacking(std::vector & blockedMatrix, - * const Eigen::MatrixXd & fullMatrix, int nrBlocks, int dimBlock) - * \param[out] blockedMatrix the result of packing fullMatrix - * \param[in] fullMatrix the matrix to be packed - * \param[in] dimBlock the dimension of the square blocks - * \param[in] nrBlocks the number of square blocks - */ -inline void symmetryPacking(std::vector & blockedMatrix, - const Eigen::MatrixXd & fullMatrix, - int dimBlock, - int nrBlocks) { - // This function packs the square block diagonal fullMatrix with nrBlocks of - // dimension dimBlock - // into a std::vector containing nrBlocks square matrices of - // dimenion dimBlock. - int j = 0; - for (int i = 0; i < nrBlocks; ++i) { - blockedMatrix.push_back(fullMatrix.block(j, j, dimBlock, dimBlock)); - j += dimBlock; - } -} - -/*! \fn inline void hermitivitize(Eigen::MatrixBase & obj_) +/*! \brief Enforce Hermitian symmetry on a matrix * \param[out] obj_ the Eigen object to be hermitivitized * \tparam Derived the numeric type of obj_ elements * @@ -246,11 +94,11 @@ inline void hermitivitize(Eigen::MatrixBase & obj_) { *rotation * matrix will be: * \f[ - * R = \begin{pmatrix} - * c_1c_3 - s_1c_2s_3 & -s_1c_3 - c_1c_2s_3 & s_2s_3 \\ - * c_1s_3 + s_1c_2c_3 & -s_1s_3 + c_1c_2c_3 & -s_2c_3 \\ - * s_1s_2 & c_1s_2 & c_2 - * \end{pmatrix} + * R = \begin{pmatrix} + * c_1c_3 - s_1c_2s_3 & -s_1c_3 - c_1c_2s_3 & s_2s_3 \\ + * c_1s_3 + s_1c_2c_3 & -s_1s_3 + c_1c_2c_3 & -s_2c_3 \\ + * s_1s_2 & c_1s_2 & c_2 + * \end{pmatrix} * \f] * Eigen's geometry module is used to calculate the rotation matrix */ @@ -265,69 +113,6 @@ inline void eulerRotation(Eigen::Matrix3d & R_, Eigen::AngleAxis(phi, Eigen::Vector3d::UnitZ()); } -/*! \brief Return value of function defined on grid at an arbitrary point - * \param[in] point where the function has to be evaluated - * \param[in] grid holds points on grid where function is known - * \param[in] function holds known function values - * - * This function finds the nearest values for the given point - * and performs a linear interpolation. - * \warning This function assumes that grid has already been sorted! - */ -inline double linearInterpolation(const double point, - const std::vector & grid, - const std::vector & function) { - // Find nearest points on grid to the arbitrary point given - size_t index = std::distance(grid.begin(), - std::lower_bound(grid.begin(), grid.end(), point)) - - 1; - - // Parameters for the interpolating line - double y_1 = function[index], y_0 = function[index - 1]; - double x_1 = grid[index], x_0 = grid[index - 1]; - double m = (y_1 - y_0) / (x_1 - x_0); - - return (m * (point - x_0) + y_0); -} - -/*! \brief Return value of function defined on grid at an arbitrary point - * \param[in] point where the function has to be evaluated - * \param[in] grid holds points on grid where function is known - * \param[in] function holds known function values - * - * This function finds the nearest values for the given point - * and performs a cubic spline interpolation. - * \warning This function assumes that grid has already been sorted! - */ -inline double splineInterpolation(const double point, - const std::vector & grid, - const std::vector & function) { - // Find nearest points on grid to the arbitrary point given - int index = - std::distance(grid.begin(), std::lower_bound(grid.begin(), grid.end(), point)); - - int imax = grid.size() - 1; - if (index <= 0) - index = 1; - if (index >= imax - 1) - index = imax - 2; - - // Parameters for the interpolating spline - Eigen::Vector4d x = (Eigen::Vector4d() << grid[index - 1], - grid[index], - grid[index + 1], - grid[index + 2]) - .finished(); - Eigen::Vector4d y = (Eigen::Vector4d() << function[index - 1], - function[index], - function[index + 1], - function[index + 2]) - .finished(); - SplineFunction s(x, y); - - return s(point); -} - /*! \brief Prints Eigen object (matrix or vector) to file * \param[in] matrix Eigen object * \param[in] fname name of the file diff --git a/src/utils/Molecule.cpp b/src/utils/Molecule.cpp index 11ae2ea6e..a6dcd2d88 100644 --- a/src/utils/Molecule.cpp +++ b/src/utils/Molecule.cpp @@ -57,7 +57,7 @@ Molecule::Molecule(int nat, atoms_(at), spheres_(sph) { rotor_ = findRotorType(); - pointGroup_ = buildGroup(0, 0, 0, 0); + pointGroup_ = Symmetry(0, 0, 0, 0); } Molecule::Molecule(int nat, @@ -75,7 +75,7 @@ Molecule::Molecule(int nat, atoms_(at), spheres_(sph) { rotor_ = findRotorType(); - pointGroup_ = buildGroup(nr_gen, gen[0], gen[1], gen[2]); + pointGroup_ = Symmetry(nr_gen, gen); } Molecule::Molecule(int nat, @@ -108,7 +108,7 @@ Molecule::Molecule(const std::vector & sph) atoms_.push_back(Atom("Dummy", "Du", charge, mass, mass, geometry_.col(i))); } rotor_ = findRotorType(); - pointGroup_ = buildGroup(0, 0, 0, 0); + pointGroup_ = Symmetry(0, 0, 0, 0); } Molecule::Molecule(const Molecule & other) { *this = other; } diff --git a/src/utils/Molecule.hpp b/src/utils/Molecule.hpp index ae3918ef8..a6016bfe2 100644 --- a/src/utils/Molecule.hpp +++ b/src/utils/Molecule.hpp @@ -47,6 +47,7 @@ class Element; namespace pcm { using utils::Atom; using utils::Sphere; +using utils::Symmetry; enum rotorType { rtAsymmetric, rtSymmetric, rtSpherical, rtLinear, rtAtom }; const std::string rotorTypeList[] = {"Asymmetric", @@ -91,7 +92,7 @@ class Molecule { */ Molecule() { rotor_ = rtAsymmetric; - pointGroup_ = buildGroup(0, 0, 0, 0); + pointGroup_ = Symmetry(0, 0, 0, 0); } /*! \brief Constructor from full molecular data * \param[in] nat number of atoms diff --git a/src/utils/RungeKutta.hpp b/src/utils/RungeKutta.hpp new file mode 100644 index 000000000..f85675609 --- /dev/null +++ b/src/utils/RungeKutta.hpp @@ -0,0 +1,197 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver 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. + * + * PCMSolver 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 PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#pragma once + +#include + +#include "Config.hpp" + +/*! \file RungeKutta.hpp */ + +namespace pcm { +namespace utils { +namespace detail { +/*! \return t1 < t2 if dt > 0 and t1 > t2 if dt < 0 with epsilon accuracy + * \note This function is part of Boost.Odeint + */ +template bool less_with_sign(T t1, T t2, T dt) { + if (dt > 0) { // return t1 < t2 + return t2 - t1 > std::numeric_limits::epsilon(); + } else { // return t1 > t2 + return t1 - t2 > std::numeric_limits::epsilon(); + } +} + +/*! \return t1 <= t2 if dt > 0 and t1 => t2 if dt < 0 with epsilon accuracy + * \note This function is part of Boost.Odeint + */ +template bool less_eq_with_sign(T t1, T t2, T dt) { + if (dt > 0) { + return t1 - t2 <= std::numeric_limits::epsilon(); + } else { + return t2 - t1 <= std::numeric_limits::epsilon(); + } +} +} // namespace detail + +/*! \brief Integrates an ODE given a stepper. + * \author Roberto Di Remigio + * \date 2018 + * \tparam Stepper integrate stepping function + * \tparam System the ODE system to integrate + * \tparam Observer function reporting integration progress + * + * The interface mimics the one of the Boost.Odeint package. + */ +template +size_t integrate_const(const Stepper & stepper, + const System & system, + typename System::StateType & f, + const double tStart, + const double tEnd, + const double dt, + const Observer & observer) { + size_t nSteps = 0; + double t = tStart; + while (detail::less_eq_with_sign(t + dt, tEnd, dt)) { + observer(f, t); + stepper.doStep(system, f, t, dt); + ++nSteps; + t = tStart + nSteps * dt; + } + observer(f, t); + return nSteps; +} + +/*! \brief 4th-order Runge-Kutta stepper. + * \author Roberto Di Remigio + * \date 2018 + * \tparam StateType the data structure used to hold the state of the ODE + * + * The interface mimics the one of the Boost.Odeint package. + */ +template struct RungeKutta4 __final { + /*! \brief Implements the "classic" 4th-order Runge-Kutta stepper. + * \author Roberto Di Remigio + * \date 2018 + * \tparam System the ODE system to integrate + * + * This function implements the "classic" Runge-Kutta stepper, with the + * following Butcher Tableau: + * + * \f[ + * \begin{array}{c|cccc} + * 0 & 0 & 0 & 0 & 0\\ + * 1/2 & 1/2 & 0 & 0 & 0\\ + * 1/2 & 0 & 1/2 & 0 & 0\\ + * 1 & 0 & 0 & 1 & 0\\ + * \hline + * & 1/6 & 1/3 & 1/3 & 1/6\\ + * \end{array} + * \f] + */ + template + void doStep(const System & system, + StateType & f, + const double t, + const double dt) const { + StateType k1, k2, k3, k4; + StateType x1, x2, x3; + + // Step 1 + system(f, k1, t); + // Step2 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x1[i] = f[i] + 0.5 * dt * k1[i]; + } + system(x1, k2, t + 0.5 * dt); + // Step 3 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x2[i] = f[i] + 0.5 * dt * k2[i]; + } + system(x2, k3, t + 0.5 * dt); + // Step 4 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x3[i] = f[i] + dt * k3[i]; + } + system(x3, k4, t + dt); + // Update + for (size_t i = 0; i < system.ODEorder(); ++i) { + f[i] += dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0; + } + } + + /*! \brief Implements the 3/8 4th-order Runge-Kutta stepper. + * \author Roberto Di Remigio + * \date 2018 + * \tparam System the ODE system to integrate + * + * This function implements the 3/8 Runge-Kutta stepper, with the + * following Butcher Tableau: + * + * \f[ + * \begin{array}{c|cccc} + * 0 & 0 & 0 & 0 & 0\\ + * 1/3 & 1/3 & 0 & 0 & 0\\ + * 2/3 & -1/3 & 1 & 0 & 0\\ + * 1 & 1 & -1 & 1 & 0\\ + * \hline + * & 1/8 & 3/8 & 3/8 & 1/8\\ + * \end{array} + * \f] + */ + template + void doStep38(const System & system, + StateType & f, + const double t, + const double dt) const { + StateType k1{}, k2{}, k3{}, k4{}; + StateType x1{}, x2{}, x3{}; + + // Step 1 + system(f, k1, t); + // Step 2 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x1[i] = f[i] + (dt / 3.0) * k1[i]; + } + system(x1, k2, t + (dt / 3.0)); + // Step 3 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x2[i] = f[i] + dt * (-1.0 / 3.0 * k1[i] + k2[i]); + } + system(x2, k3, t + dt * (2.0 / 3.0)); + // Step 4 + for (size_t i = 0; i < system.ODEorder(); ++i) { + x3[i] = f[i] + dt * (k1[i] - k2[i] + k3[i]); + } + system(x3, k4, t + dt); + // Update + for (size_t i = 0; i < system.ODEorder(); ++i) { + f[i] += dt * ((1.0 / 8.0) * k1[i] + (3.0 / 8.0) * k2[i] + (3.0 / 8.0) * k3[i] + + (1.0 / 8.0) * k4[i]); + } + } +}; +} // namespace utils +} // namespace pcm diff --git a/src/utils/RungeKutta4.cpp b/src/utils/RungeKutta4.cpp deleted file mode 100644 index c35a40502..000000000 --- a/src/utils/RungeKutta4.cpp +++ /dev/null @@ -1,134 +0,0 @@ -/* - * The Runge-Kutta source code is Copyright(c) 2010 John Burkardt. - * - * This Runge-Kutta ODE solver 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. - * - * This Runge-Kutts ODE Solver 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 the code. If not, see . - * - * See also: - * https://people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html - */ - -#include "rk4.hpp" - -#include -#include -#include -#include - -namespace pcm { -namespace utils { -using namespace std; - -double rk4(double t0, double u0, double dt, double f(double t, double u)) { - double f0; - double f1; - double f2; - double f3; - double t1; - double t2; - double t3; - double u; - double u1; - double u2; - double u3; - // - // Get four sample values of the derivative. - // - f0 = f(t0, u0); - - t1 = t0 + dt / 2.0; - u1 = u0 + dt * f0 / 2.0; - f1 = f(t1, u1); - - t2 = t0 + dt / 2.0; - u2 = u0 + dt * f1 / 2.0; - f2 = f(t2, u2); - - t3 = t0 + dt; - u3 = u0 + dt * f2; - f3 = f(t3, u3); - // - // Combine to estimate the solution at time T0 + DT. - // - u = u0 + dt * (f0 + 2.0 * f1 + 2.0 * f2 + f3) / 6.0; - - return u; -} - -double * rk4vec(double t0, - int m, - double u0[], - double dt, - double * f(double t, int m, double u[])) - -{ - double * f0; - double * f1; - double * f2; - double * f3; - int i; - double t1; - double t2; - double t3; - double * u; - double * u1; - double * u2; - double * u3; - // - // Get four sample values of the derivative. - // - f0 = f(t0, m, u0); - - t1 = t0 + dt / 2.0; - u1 = new double[m]; - for (i = 0; i < m; i++) { - u1[i] = u0[i] + dt * f0[i] / 2.0; - } - f1 = f(t1, m, u1); - - t2 = t0 + dt / 2.0; - u2 = new double[m]; - for (i = 0; i < m; i++) { - u2[i] = u0[i] + dt * f1[i] / 2.0; - } - f2 = f(t2, m, u2); - - t3 = t0 + dt; - u3 = new double[m]; - for (i = 0; i < m; i++) { - u3[i] = u0[i] + dt * f2[i]; - } - f3 = f(t3, m, u3); - // - // Combine them to estimate the solution. - // - u = new double[m]; - for (i = 0; i < m; i++) { - u[i] = u0[i] + dt * (f0[i] + 2.0 * f1[i] + 2.0 * f2[i] + f3[i]) / 6.0; - } - // - // Free memory. - // - delete[] f0; - delete[] f1; - delete[] f2; - delete[] f3; - delete[] u1; - delete[] u2; - delete[] u3; - - return u; -} -} // namespace utils -} // namespace pcm diff --git a/src/utils/RungeKutta4.hpp b/src/utils/RungeKutta4.hpp deleted file mode 100644 index 8759168a2..000000000 --- a/src/utils/RungeKutta4.hpp +++ /dev/null @@ -1,82 +0,0 @@ -/* - * The Runge-Kutta source code is Copyright(c) 2013 John Burkardt. - * - * This Runge-Kutta ODE solver 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. - * - * This Runge-Kutts ODE Solver 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 the code. If not, see . - * - * See also: - * https://people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html - */ - -/*! \file RungeKutta4.hpp */ - -namespace pcm { -namespace utils { -/*! Right-hand side for a scalar ordinary-differential equation */ -typedef std::function ScalarODERHS; - -/*! \brief RK4 takes one Runge-Kutta step for a scalar ODE. - * \author John Burkardt, Roberto Di Remigio - * \param[in] t0 the current time. - * \param[in] u0 the solution estimate at the current time. - * \param[in] dt the time step. - * \param[in] F a function which evaluates the derivative, or right hand side of the - * problem. - * \return the fourth-order Runge-Kutta solution estimate at time t0+dt. - * - * It is assumed that an initial value problem, of the form - * - * du/dt = f( t, u ) - * u(t0) = u0 - * - * is being solved. - * - * If the user can supply current values of t, u, a stepsize dt, and a - * function to evaluate the derivative, this function can compute the - * fourth-order Runge Kutta estimate to the solution at time t+dt. - */ -double rk4(double t0, double u0, double dt, double f(double t, double u)); - -/*! Right-hand side for a vector ordinary-differential equation */ -typedef std::function VectorODERHS; - -/*! - * \brief Takes one Runge-Kutta step for a vector ODE. - * \author John Burkardt, Roberto Di Remigio - * \param[in] t0 the current time. - * \param[in] n the spatial dimension. - * \param[in] u0[M] the solution estimate at the current time. - * \param[in] dt the time step. - * \param[in] F a function which evaluates the derivative, or right hand side of the - * problem. - * \return the fourth-order Runge-Kutta solution estimate at time t0+dt. - * - * It is assumed that an initial value problem, of the form - * - * du/dt = f ( t, u ) - * u(t0) = u0 - * - * is being solved. - * - * If the user can supply current values of t, u, a stepsize dt, and a - * function to evaluate the derivative, this function can compute the - * fourth-order Runge Kutta estimate to the solution at time t+dt. - */ -double * rk4vec(double t0, - int n, - double u0[], - double dt, - double * f(double t, int n, double u[])); -} // namespace utils -} // namespace pcm diff --git a/src/utils/SplineFunction.hpp b/src/utils/SplineFunction.hpp index 75b73008b..6030b88a4 100644 --- a/src/utils/SplineFunction.hpp +++ b/src/utils/SplineFunction.hpp @@ -24,21 +24,24 @@ #pragma once #include +#include #include "Config.hpp" #include #include -/*! \file SplineFunction.hpp - * \class SplineFunction +/*! \file SplineFunction.hpp */ + +namespace pcm { +namespace utils { +/*! \class SplineFunction * \brief Spline interpolation of a function * \author Roberto Di Remigio * \date 2015 * * Taken from StackOverflow http://stackoverflow.com/a/29825204/2528668 */ - class SplineFunction __final { private: typedef Eigen::Spline CubicSpline; @@ -94,3 +97,43 @@ class SplineFunction __final { } #endif /*HAS_CXX11_LAMBDA */ }; + +/*! \brief Return value of function defined on grid at an arbitrary point + * \param[in] point where the function has to be evaluated + * \param[in] grid holds points on grid where function is known + * \param[in] function holds known function values + * + * This function finds the nearest values for the given point + * and performs a cubic spline interpolation. + * \warning This function assumes that grid has already been sorted! + */ +inline double splineInterpolation(const double point, + const std::vector & grid, + const std::vector & function) { + // Find nearest points on grid to the arbitrary point given + int index = + std::distance(grid.begin(), std::lower_bound(grid.begin(), grid.end(), point)); + + int imax = grid.size() - 1; + if (index <= 0) + index = 1; + if (index >= imax - 1) + index = imax - 2; + + // Parameters for the interpolating spline + Eigen::Vector4d x = (Eigen::Vector4d() << grid[index - 1], + grid[index], + grid[index + 1], + grid[index + 2]) + .finished(); + Eigen::Vector4d y = (Eigen::Vector4d() << function[index - 1], + function[index], + function[index + 1], + function[index + 2]) + .finished(); + SplineFunction s(x, y); + + return s(point); +} +} // namespace utils +} // namespace pcm diff --git a/src/utils/Symmetry.cpp b/src/utils/Symmetry.cpp index 0032d5d9a..4cbc61fa5 100644 --- a/src/utils/Symmetry.cpp +++ b/src/utils/Symmetry.cpp @@ -23,27 +23,93 @@ #include "Symmetry.hpp" +#include + #include "Config.hpp" #include -/* Indexing of symmetry operations and their mapping to a bitstring: - * zyx Parity - * 0 000 E 1.0 - * 1 001 Oyz -1.0 - * 2 010 Oxz -1.0 - * 3 011 C2z 1.0 - * 4 100 Oxy -1.0 - * 5 101 C2y 1.0 - * 6 110 C2x 1.0 - * 7 111 i -1.0 - */ +#include "MathUtils.hpp" -Symmetry buildGroup(int _nr_gen, int _gen1, int _gen2, int _gen3) { - int gen[3]; - gen[0] = _gen1; - gen[1] = _gen2; - gen[2] = _gen3; +namespace pcm { +namespace utils { +void symmetryBlocking(Eigen::MatrixXd & matrix, + PCMSolverIndex cavitySize, + PCMSolverIndex ntsirr, + int nr_irrep) { + // u is the character table for the group (t in the paper) + Eigen::MatrixXd u = Eigen::MatrixXd::Zero(nr_irrep, nr_irrep); + for (int i = 0; i < nr_irrep; ++i) { + for (int j = 0; j < nr_irrep; ++j) { + u(i, j) = parity(i & j); + } + } + // Naming of indices: + // a, b, c, d run over the total size of the cavity (N) + // i, j, k, l run over the number of irreps (n) + // p, q, r, s run over the irreducible size of the cavity (N/n) + // Instead of forming U (T in the paper) and then perform the dense + // multiplication, we multiply block-by-block using just the u matrix. + // matrix = U * matrix * Ut; U * Ut = Ut * U = id + // First half-transformation, i.e. first_half = matrix * Ut + Eigen::MatrixXd first_half = Eigen::MatrixXd::Zero(cavitySize, cavitySize); + for (int i = 0; i < nr_irrep; ++i) { + int ioff = i * ntsirr; + for (int k = 0; k < nr_irrep; ++k) { + int koff = k * ntsirr; + for (int j = 0; j < nr_irrep; ++j) { + int joff = j * ntsirr; + double ujk = u(j, k) / nr_irrep; + for (int p = 0; p < ntsirr; ++p) { + int a = ioff + p; + for (int q = 0; q < ntsirr; ++q) { + int b = joff + q; + int c = koff + q; + first_half(a, c) += matrix(a, b) * ujk; + } + } + } + } + } + // Second half-transformation, i.e. matrix = U * first_half + matrix.setZero(cavitySize, cavitySize); + for (int i = 0; i < nr_irrep; ++i) { + int ioff = i * ntsirr; + for (int k = 0; k < nr_irrep; ++k) { + int koff = k * ntsirr; + for (int j = 0; j < nr_irrep; ++j) { + int joff = j * ntsirr; + double uij = u(i, j); + for (int p = 0; p < ntsirr; ++p) { + int a = ioff + p; + int b = joff + p; + for (int q = 0; q < ntsirr; ++q) { + int c = koff + q; + matrix(a, c) += uij * first_half(b, c); + } + } + } + } + } + // Traverse the matrix and discard numerical zeros + for (PCMSolverIndex a = 0; a < cavitySize; ++a) { + for (PCMSolverIndex b = 0; b < cavitySize; ++b) { + if (numericalZero(matrix(a, b))) { + matrix(a, b) = 0.0; + } + } + } +} - return Symmetry(_nr_gen, gen); +void symmetryPacking(std::vector & blockedMatrix, + const Eigen::MatrixXd & fullMatrix, + int dimBlock, + int nrBlocks) { + int j = 0; + for (int i = 0; i < nrBlocks; ++i) { + blockedMatrix.push_back(fullMatrix.block(j, j, dimBlock, dimBlock)); + j += dimBlock; + } } +} // namespace utils +} // namespace pcm diff --git a/src/utils/Symmetry.hpp b/src/utils/Symmetry.hpp index db19b05f3..a7676ab3c 100644 --- a/src/utils/Symmetry.hpp +++ b/src/utils/Symmetry.hpp @@ -24,59 +24,151 @@ #pragma once #include +#include #include +#include #include "Config.hpp" -/*! \file Symmetry.hpp - * \class Symmetry +#include + +/*! \file Symmetry.hpp */ + +namespace pcm { +namespace utils { +/*! \class Symmetry * \brief Contains very basic info about symmetry (only Abelian groups) * \author Roberto Di Remigio - * \date 2014 + * \date 2014, 2018 + * \note C1 is built as Symmetry C1(0, 0, 0, 0); * - * Just a wrapper around a vector containing the generators of the group + * Indexing of symmetry operations and their mapping to a bitstring: + * zyx Parity + * 0 000 E 1.0 + * 1 001 Oyz -1.0 + * 2 010 Oxz -1.0 + * 3 011 C2z 1.0 + * 4 100 Oxy -1.0 + * 5 101 C2y 1.0 + * 6 110 C2x 1.0 + * 7 111 i -1.0 */ - -class Symmetry { -private: - /*! - * Number of generators - */ - int nrGenerators_; - /*! - * Generators - */ - int generators_[3]; - /*! - * Number of irreps - */ - int nrIrrep_; - +class Symmetry __final { public: - /*! \brief Default constructor sets up C1 point group - */ - Symmetry() : nrGenerators_(0) { +#ifdef HAS_CXX11 + Symmetry() + : nrGenerators_(0), + nrIrrep_(static_cast(std::pow(2.0, nrGenerators_))), + generators_({{0, 0, 0}}) {} + Symmetry(int nr_gen, int gen1, int gen2, int gen3) + : nrGenerators_(nr_gen), + nrIrrep_(static_cast(std::pow(2.0, nrGenerators_))), + generators_({{gen1, gen2, gen3}}) {} +#else + Symmetry() + : nrGenerators_(0), nrIrrep_(static_cast(std::pow(2.0, nrGenerators_))) { std::fill(generators_, generators_ + 3, 0); - nrIrrep_ = int(std::pow(2.0, nrGenerators_)); } - Symmetry(int nr_gen, int gen[3]) : nrGenerators_(nr_gen) { - // Transfer the passed generators array into generators_ - std::copy(gen, gen + nrGenerators_, generators_); - // We can now initialize the number of irreps - nrIrrep_ = int(std::pow(2.0, nrGenerators_)); + Symmetry(int nr_gen, int gen1, int gen2, int gen3) + : nrGenerators_(nr_gen), + nrIrrep_(static_cast(std::pow(2.0, nrGenerators_))), + { + generators_[0] = gen1; + generators_[1] = gen2; + generators_[2] = gen3; } - Symmetry(const Symmetry & other) - : nrGenerators_(other.nrGenerators_), nrIrrep_(other.nrIrrep_) { - std::copy(other.generators_, other.generators_ + nrGenerators_, generators_); +#endif + Symmetry(int nr_gen, int gen[3]) + : nrGenerators_(nr_gen), + nrIrrep_(static_cast(std::pow(2.0, nrGenerators_))) { + generators_[0] = gen[0]; + generators_[1] = gen[1]; + generators_[2] = gen[2]; } - ~Symmetry() {} + int nrGenerators() const { return nrGenerators_; } - int generators(int i) const { return generators_[i]; } int nrIrrep() const { return nrIrrep_; } + int generators(size_t i) const { return generators_[i]; } + +private: + /*! Number of generators */ + int nrGenerators_; + /*! Number of irreps */ + int nrIrrep_; + /*! Generators */ + pcm::array generators_; }; -/*! Builds Symmetry object. +/*! \brief Transform matrix to block diagonal form by symmetry + * \param[out] matrix the matrix to be block-diagonalized + * \param[in] cavitySize the size of the cavity (size of the matrix) + * \param[in] ntsirr the size of the irreducible portion of the cavity (size of + * the blocks) + * \param[in] nr_irrep the number of irreducible representations (number of + * blocks) + * + * This function implements the symmetry-blocking of the PCM matrix due to + * point group symmetry as reported in: + * L. Frediani, R. Cammi, C. S. Pomelli, J. Tomasi and K. Ruud, J. Comput.Chem. 25, + * 375 (2003) + */ +void symmetryBlocking(Eigen::MatrixXd & matrix, + PCMSolverIndex cavitySize, + PCMSolverIndex ntsirr, + int nr_irrep); + +/*! \brief Packs symmetry blocked matrix, i.e. only stores non-zero blocks + * \param[out] blockedMatrix the result of packing fullMatrix + * \param[in] fullMatrix the matrix to be packed + * \param[in] dimBlock the dimension of the square blocks + * \param[in] nrBlocks the number of square blocks + * + * This function packs the square block diagonal fullMatrix with nrBlocks of + * dimension dimBlock into a std::vector containing nrBlocks + * square matrices of dimenion dimBlock. + */ +void symmetryPacking(std::vector & blockedMatrix, + const Eigen::MatrixXd & fullMatrix, + int dimBlock, + int nrBlocks); + +/*! \brief Calculate the parity of the bitset + * \param[in] bitrep a bitset + * \tparam nBits length of the input bitset + * + * Calculate the parity of the bitset as defined by: + * bitrep[0] XOR bitrep[1] XOR ... XOR bitrep[nBits-1] + */ +template inline int parity(std::bitset bitrep) { + int parity = 0; + for (size_t i = 0; i < bitrep.size(); ++i) { + parity ^= bitrep[i]; + } + return parity; +} + +/*! \brief Returns parity of input integer. + * \param[in] i an integer, usually an index for an irrep or a symmetry operation + * + * The parity is defined as the result of using XOR on the bitrep + * of the given integer. For example: + * 2 -> 010 -> 0^1^0 = 1 -> -1.0 + * 6 -> 110 -> 1^1^0 = 0 -> 1.0 * - * \note C1 is built as Symmetry C1 = buildGroup(0, 0, 0, 0); + * It can also be interpreted as the action of a given operation on the + * Cartesian axes: + * zyx Parity + * 0 000 E 1.0 + * 1 001 Oyz -1.0 + * 2 010 Oxz -1.0 + * 3 011 C2z 1.0 + * 4 100 Oxy -1.0 + * 5 101 C2y 1.0 + * 6 110 C2x 1.0 + * 7 111 i -1.0 */ -Symmetry buildGroup(int _nr_gen, int _gen1, int _gen2, int _gen3); +inline double parity(unsigned int i) { + return (parity(std::bitset<3>(i)) ? -1.0 : 1.0); +} +} // namespace utils +} // namespace pcm diff --git a/tests/TestingMolecules.hpp b/tests/TestingMolecules.hpp index ce51f7d3e..e7e1e8a32 100644 --- a/tests/TestingMolecules.hpp +++ b/tests/TestingMolecules.hpp @@ -82,38 +82,38 @@ template Molecule dummy(double radius, const Eigen::Vector3d & cente Symmetry pGroup; switch (group) { case (pgC1): - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; case (pgC2): // C2 as generated by C2z - pGroup = buildGroup(1, 3, 0, 0); + pGroup = Symmetry(1, 3, 0, 0); break; case (pgCs): // Cs as generated by Oyz - pGroup = buildGroup(1, 1, 0, 0); + pGroup = Symmetry(1, 1, 0, 0); break; case (pgCi): // Ci as generated by i - pGroup = buildGroup(1, 7, 0, 0); + pGroup = Symmetry(1, 7, 0, 0); break; case (pgD2): // D2 as generated by C2z and C2x - pGroup = buildGroup(2, 3, 6, 0); + pGroup = Symmetry(2, 3, 6, 0); break; case (pgC2v): // C2v as generated by Oyz and Oxz - pGroup = buildGroup(2, 1, 2, 0); + pGroup = Symmetry(2, 1, 2, 0); break; case (pgC2h): // C2h as generated by Oxy and i - pGroup = buildGroup(2, 4, 7, 0); + pGroup = Symmetry(2, 4, 7, 0); break; case (pgD2h): // D2h as generated by Oxy, Oxz and Oyz - pGroup = buildGroup(3, 4, 2, 1); + pGroup = Symmetry(3, 4, 2, 1); break; default: - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; } @@ -159,7 +159,7 @@ Molecule NH3() { spheres.push_back(sph4); // C1 - Symmetry pGroup = buildGroup(0, 0, 0, 0); + Symmetry pGroup(0, 0, 0, 0); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; @@ -197,14 +197,14 @@ template Molecule H3() { Symmetry pGroup; switch (group) { case (pgC1): - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; case (pgC2v): // C2v as generated by Oyz and Oxz - pGroup = buildGroup(2, 1, 2, 0); + pGroup = Symmetry(2, 1, 2, 0); break; default: - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; } @@ -235,7 +235,7 @@ Molecule H2() { spheres.push_back(sph2); spheres.push_back(sph3); - Symmetry pGroup = buildGroup(0, 0, 0, 0); + Symmetry pGroup(0, 0, 0, 0); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; @@ -274,38 +274,38 @@ template Molecule CO2() { Symmetry pGroup; switch (group) { case (pgC1): - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; case (pgC2): // C2 as generated by C2z - pGroup = buildGroup(1, 3, 0, 0); + pGroup = Symmetry(1, 3, 0, 0); break; case (pgCs): // Cs as generated by Oyz - pGroup = buildGroup(1, 1, 0, 0); + pGroup = Symmetry(1, 1, 0, 0); break; case (pgCi): // Ci as generated by i - pGroup = buildGroup(1, 7, 0, 0); + pGroup = Symmetry(1, 7, 0, 0); break; case (pgD2): // D2 as generated by C2z and C2x - pGroup = buildGroup(2, 3, 6, 0); + pGroup = Symmetry(2, 3, 6, 0); break; case (pgC2v): // C2v as generated by Oyz and Oxz - pGroup = buildGroup(2, 1, 2, 0); + pGroup = Symmetry(2, 1, 2, 0); break; case (pgC2h): // C2h as generated by Oxy and i - pGroup = buildGroup(2, 4, 7, 0); + pGroup = Symmetry(2, 4, 7, 0); break; case (pgD2h): // D2h as generated by Oxy, Oxz and Oyz - pGroup = buildGroup(3, 4, 2, 1); + pGroup = Symmetry(3, 4, 2, 1); break; default: - pGroup = buildGroup(0, 0, 0, 0); + pGroup = Symmetry(0, 0, 0, 0); break; } @@ -348,7 +348,7 @@ Molecule CH3() { spheres.push_back(sph4); // Cs as generated by Oxy - Symmetry pGroup = buildGroup(1, 4, 0, 0); + Symmetry pGroup(1, 4, 0, 0); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; @@ -399,7 +399,7 @@ Molecule C2H4() { spheres.push_back(sph6); // D2h as generated by Oxy, Oxz, Oyz - Symmetry pGroup = buildGroup(3, 4, 2, 1); + Symmetry pGroup(3, 4, 2, 1); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; @@ -496,7 +496,7 @@ Molecule C6H6() { spheres.push_back(sph11); spheres.push_back(sph12); - Symmetry pGroup = buildGroup(0, 0, 0, 0); + Symmetry pGroup(0, 0, 0, 0); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; @@ -533,7 +533,7 @@ Molecule H2O() { spheres.push_back(sph2); spheres.push_back(sph3); - Symmetry pGroup = buildGroup(0, 0, 0, 0); + Symmetry pGroup(0, 0, 0, 0); return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); }; diff --git a/tests/green/CMakeLists.txt b/tests/green/CMakeLists.txt index e280458ea..58236b785 100644 --- a/tests/green/CMakeLists.txt +++ b/tests/green/CMakeLists.txt @@ -57,13 +57,13 @@ add_Catch_test( green_anisotropic_liquid ) -# green_sharp_diffuse.cpp test +# green_spherical_sharp.cpp test add_Catch_test( NAME - green_sharp_diffuse + green_spherical_sharp LABELS green - green_sharp_diffuse + green_spherical_sharp COST 100.0 ) From 8b2f3a9185e4830d1a07d7479f8c1e5b243c9e49 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Tue, 27 Mar 2018 10:16:37 -0400 Subject: [PATCH 3/3] Updtate reference numbers in test --- tests/green/green_spherical_diffuse.cpp | 70 ++++++++++++------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/tests/green/green_spherical_diffuse.cpp b/tests/green/green_spherical_diffuse.cpp index f634ccc5a..7d018bceb 100644 --- a/tests/green/green_spherical_diffuse.cpp +++ b/tests/green/green_spherical_diffuse.cpp @@ -68,7 +68,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero(); SphericalDiffuse<> gf(eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0204749331147813275; + double value = 0.0204749331147992819; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -79,7 +79,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.488060362514614432; + double value = 0.487473605103865837; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -91,7 +91,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.0124999765635167015; + double derProbe = -0.0124999765637075211; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -103,7 +103,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.290662079075743041; + double derProbe = -0.290662955233778053; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -115,7 +115,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the source " "point inside the droplet is") { - double derSource = 0.00937504723103055326; + double derSource = 0.00937504723101320603; double gf_derSource = gf.derivativeSource(sourceNormal1, source1, probe1); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -127,7 +127,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = -0.938830970128590181; + double derSource = -0.901954472843635724; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -144,7 +144,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ (Eigen::Vector3d() << 25.0, 0.0, 0.0).finished(); SphericalDiffuse<> gf(eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0173845118865096904; + double value = 0.0173845118865097424; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -155,7 +155,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.501844654162973303; + double value = 0.501385603574595162; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -167,7 +167,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.012483438527558649; + double derProbe = -0.0124834385275933435; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -179,7 +179,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.291257922496734878; + double derProbe = -0.291262383059764929; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -191,7 +191,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the source " "point inside the droplet is") { - double derSource = 0.0124835646475064677; + double derSource = 0.012483564647454426; double gf_derSource = gf.derivativeSource(sourceNormal1, source1, probe1); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -203,7 +203,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = 5.01504540067754245; + double derSource = 2.720604738008503; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -221,7 +221,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ SphericalDiffuse gf( eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0204265162808782985; + double value = 0.0204265162808963327; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -232,7 +232,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.428395676083591359; + double value = 0.42831756921770181; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -244,7 +244,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.0124999769449823939; + double derProbe = -0.0124999769448609632; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -256,7 +256,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.288385047571004804; + double derProbe = -0.288392889131738883; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -280,7 +280,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = 3.22957321418654297; + double derSource = -1.51576208806569745; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -298,7 +298,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ SphericalDiffuse gf( eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0173358023608628994; + double value = 0.017335802360862948; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -309,7 +309,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.471193749836353093; + double value = 0.470725113611886292; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -321,7 +321,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.0124834504454732209; + double derProbe = -0.0124834504460630269; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -333,7 +333,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.290867341330436346; + double derProbe = -0.290870450536662162; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -345,7 +345,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the source " "point inside the droplet is") { - double derSource = 0.0124835522706361057; + double derSource = 0.012483552270584064; double gf_derSource = gf.derivativeSource(sourceNormal1, source1, probe1); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -357,7 +357,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = 4.9830666203676266; + double derSource = 3.75908425086779463; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -376,7 +376,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ SphericalDiffuse gf( eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0204465601345453427; + double value = 0.0204465601345278428; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -387,7 +387,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.546459075558889285; + double value = 0.546440841068159378; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -399,7 +399,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.0124999769317464537; + double derProbe = -0.0124999769314168563; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -411,7 +411,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.292047321249766512; + double derProbe = -0.292045737975143993; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -423,7 +423,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the source " "point inside the droplet is") { - double derSource = 0.00937504649521289646; + double derSource = 0.00937504649523024369; double gf_derSource = gf.derivativeSource(sourceNormal1, source1, probe1); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -435,7 +435,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = -2.1857233870498094; + double derSource = -1.60173403890240262; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -453,7 +453,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ SphericalDiffuse gf( eps1, eps2, width, sphereRadius, sphereCenter, maxL); THEN("the value of the Green's function inside the droplet is") { - double value = 0.0173558529449408631; + double value = 0.0173558563009676217; double gf_value = gf.kernelS(source1, probe1); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -464,7 +464,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ REQUIRE(value == Approx(gf_value)); } AND_THEN("the value of the Green's function outside the droplet is") { - double value = 0.372262125530776311; + double value = 0.37238664692877399; double gf_value = gf.kernelS(source2, probe2); INFO("ref_value = " << std::setprecision( std::numeric_limits::digits10) @@ -476,7 +476,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the probe " "point inside the droplet is") { - double derProbe = -0.0124834509424368023; + double derProbe = -0.0124834500456541542; double gf_derProbe = gf.derivativeProbe(probeNormal1, source1, probe1); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -488,7 +488,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "probe point outside the droplet is") { - double derProbe = -0.27916006373779334; + double derProbe = -0.279220523679291066; double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2); INFO("ref_derProbe = " << std::setprecision(std::numeric_limits::digits10) @@ -500,7 +500,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } THEN("the value of the Green's function directional derivative wrt the source " "point inside the droplet is") { - double derSource = 0.0124835526376168571; + double derSource = 0.0124835526869870872; double gf_derSource = gf.derivativeSource(sourceNormal1, source1, probe1); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10) @@ -512,7 +512,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ } AND_THEN("the value of the Green's function directional derivative wrt the " "source point outside the droplet is") { - double derSource = 3.4616694327527231; + double derSource = 4.79631817304426722; double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2); INFO("ref_derSource = " << std::setprecision(std::numeric_limits::digits10)