diff --git a/CHANGELOG.md b/CHANGELOG.md
index efbeae3..23dabbc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,8 +2,6 @@
## [Unreleased]
-## [Version 2.0.0-alpha1] - 2019-04-XY
-
We have introduced a number of breaking changes, motivated by the need to
modernize the library.
@@ -14,6 +12,10 @@ modernize the library.
- **BREAKING** We require a fully compliant C++14 compiler.
- **BREAKING** We require CMake 3.9
- Project dependencies are now managed as a superbuild.
+- 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.
### Removed
diff --git a/src/green/InterfacesImpl.hpp b/src/green/InterfacesImpl.hpp
new file mode 100644
index 0000000..94b614d
--- /dev/null
+++ b/src/green/InterfacesImpl.hpp
@@ -0,0 +1,277 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 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
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#include "utils/MathUtils.hpp"
+#include "utils/RungeKutta.hpp"
+#include "utils/SplineFunction.hpp"
+
+/*! \file InterfacesImpl.hpp */
+
+namespace pcm {
+namespace green {
+namespace detail {
+/*! \brief Abstract class for a 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 std::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 std::function(const double)> ProfileEvaluator;
+
+/*! \class LnTransformedRadial
+ * \brief system of ln-transformed first-order radial differential equations
+ * \author Roberto Di Remigio
+ * \date 2015
+ *
+ * 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::green::detail::ODESystem<2> {
+public:
+ /*! Type of the state vector of the ODE */
+ 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 */
+ const ProfileEvaluator eval_;
+ /*! Angular momentum */
+ 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
+ double eps = 0.0, epsPrime = 0.0;
+ std::tie(eps, epsPrime) = eval_(std::exp(y));
+ if (utils::numericalZero(eps))
+ PCMSOLVER_ERROR("Division by zero!");
+ double gamma_epsilon = std::exp(y) * epsPrime / eps;
+ // System of equations is defined here
+ drhodr[0] = rho[1];
+ drhodr[1] = -rho[1] * (rho[1] + 1.0 + gamma_epsilon) + l_ * (l_ + 1);
+ }
+};
+} // namespace detail
+
+using detail::ProfileEvaluator;
+
+/*! \class RadialFunction
+ * \brief represents solutions to the radial 2nd order ODE
+ * \author Roberto Di Remigio
+ * \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 RadialFunction final {
+public:
+ RadialFunction() : L_(0), y_min_(0.0), y_max_(0.0), y_sign_(1) {}
+ RadialFunction(int l,
+ double ymin,
+ double ymax,
+ double ystep,
+ const ProfileEvaluator & eval)
+ : L_(l),
+ y_min_(ymin),
+ y_max_(ymax),
+ y_sign_(pcm::utils::sign(y_max_ - y_min_)) {
+ compute(ystep, eval);
+ }
+ std::tuple operator()(double point) const {
+ return std::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) {
+ // 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;
+ }
+
+private:
+ typedef typename ODE::StateType StateType;
+ /*! Angular momentum of the solution */
+ int L_;
+ /*! 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 */
+ std::array, 3> function_;
+ /*! Reports progress of differential equation integrator */
+ 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 radial solution
+ * \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 (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_.
+ */
+ size_t compute(const double step, const ProfileEvaluator & eval) {
+ ODE system(eval, L_);
+ // 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_;
+ } else { // omega-type solution
+ init[0] = y_sign_ * (L_ + 1) * y_min_;
+ init[1] = y_sign_ * (L_ + 1);
+ }
+ pcm::utils::RungeKutta4 stepper;
+ size_t nSteps =
+ pcm::utils::integrate_const(stepper,
+ system,
+ init,
+ y_min_,
+ y_max_,
+ y_sign_ * step,
+ std::bind(&RadialFunction::push_back,
+ this,
+ std::placeholders::_1,
+ std::placeholders::_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) {
+ for (auto & comp : function_) {
+ std::reverse(comp.begin(), comp.end());
+ }
+ // clang-format on
+ }
+ return nSteps;
+ }
+ /*! \brief Returns value of function at given point
+ * \param[in] point evaluation point
+ *
+ * We first check if point is below y_min_, if yes we use
+ * the asymptotic form.
+ */
+ double function_impl(double point) const {
+ 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 {
+ val = utils::splineInterpolation(point, function_[0], function_[1]);
+ }
+ return val;
+ }
+ /*! \brief Returns value of 1st derivative of function at given point
+ * \param[in] point evaluation point
+ *
+ * Below y_min_, the asymptotic form is used. Otherwise we interpolate.
+ */
+ double derivative_impl(double point) const {
+ 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 {
+ val = utils::splineInterpolation(point, function_[0], function_[2]);
+ }
+ return val;
+ }
+};
+
+/*! \brief Write contents of a RadialFunction to file
+ * \param[in] f solution to the radial 2nd order ODE to print
+ * \param[in] fname name of the file
+ * \author Roberto Di Remigio
+ * \date 2018
+ * \tparam ODE system of 1st order ODEs replacing the 2nd order ODE
+ */
+template
+void writeToFile(RadialFunction & f, const std::string & fname) {
+ std::ofstream fout;
+ fout.open(fname.c_str());
+ fout << f << std::endl;
+ fout.close();
+}
+} // namespace green
+} // namespace pcm
diff --git a/src/green/SphericalDiffuse.cpp b/src/green/SphericalDiffuse.cpp
new file mode 100644
index 0000000..805fa99
--- /dev/null
+++ b/src/green/SphericalDiffuse.cpp
@@ -0,0 +1,404 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 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:
+ */
+
+#include "SphericalDiffuse.hpp"
+
+#include
+
+#include
+
+// Has to be included here
+#include "InterfacesImpl.hpp"
+
+#include "GreenData.hpp"
+#include "GreensFunction.hpp"
+#include "cavity/Element.hpp"
+#include "utils/Legendre.hpp"
+#include "utils/MathUtils.hpp"
+
+namespace pcm {
+namespace green {
+template
+double SphericalDiffuse::coefficientCoulomb(
+ const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const {
+ // Obtain coefficient for the separation of the Coulomb singularity
+ return this->coefficient_impl(source, probe);
+}
+
+template
+double SphericalDiffuse::Coulomb(
+ const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const {
+ double r12 = (source - probe).norm();
+
+ // Obtain coefficient for the separation of the Coulomb singularity
+ return (1.0 / (this->coefficient_impl(source, probe) * r12));
+}
+
+template
+double SphericalDiffuse::imagePotential(
+ const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const {
+ // Obtain coefficient for the separation of the Coulomb singularity
+ double Cr12 = this->coefficient_impl(source, probe);
+
+ double gr12 = 0.0;
+ for (int L = 1; L <= maxLGreen_; ++L) {
+ gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
+ }
+
+ return gr12;
+}
+
+template
+double SphericalDiffuse::coefficientCoulombDerivative(
+ const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const {
+ return threePointStencil(
+ std::bind(&SphericalDiffuse::coefficientCoulomb,
+ this,
+ std::placeholders::_1,
+ std::placeholders::_2),
+ p2,
+ p1,
+ direction,
+ this->delta_);
+}
+
+template
+double SphericalDiffuse::CoulombDerivative(
+ const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const {
+ return threePointStencil(std::bind(&SphericalDiffuse::Coulomb,
+ this,
+ std::placeholders::_1,
+ std::placeholders::_2),
+ p2,
+ p1,
+ direction,
+ this->delta_);
+}
+
+template
+double SphericalDiffuse::imagePotentialDerivative(
+ const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const {
+ return threePointStencil(
+ std::bind(&SphericalDiffuse::imagePotential,
+ this,
+ std::placeholders::_1,
+ std::placeholders::_2),
+ p2,
+ p1,
+ direction,
+ this->delta_);
+}
+
+template
+std::tuple SphericalDiffuse::epsilon(
+ const Eigen::Vector3d & point) const {
+ return this->profile_((point + this->origin_).norm());
+}
+
+template
+void SphericalDiffuse::toFile(const std::string & prefix) {
+ std::string tmp;
+ prefix.empty() ? tmp = prefix : tmp = prefix + "-";
+ writeToFile(zetaC_, tmp + "zetaC.dat");
+ writeToFile(omegaC_, tmp + "omegaC.dat");
+ for (int L = 1; L <= maxLGreen_; ++L) {
+ 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");
+ }
+}
+
+template
+Stencil SphericalDiffuse::operator()(Stencil * sp,
+ Stencil * pp) const {
+ // Transfer raw arrays to Eigen vectors using the Map type
+ Eigen::Map> source(sp), probe(pp);
+
+ // Obtain coefficient for the separation of the Coulomb singularity
+ double Cr12 = this->coefficient_impl(source, probe);
+
+ double gr12 = 0.0;
+ for (int L = 0; L <= maxLGreen_; ++L) {
+ gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
+ }
+ double r12 = (source - probe).norm();
+
+ double gf = (1.0 / (Cr12 * r12) + gr12);
+ return gf;
+}
+
+template
+double SphericalDiffuse::kernelD_impl(
+ const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const {
+ double eps_r2 = 0.0;
+ // Shift p2 by origin_
+ std::tie(eps_r2, std::ignore) = this->epsilon(p2);
+
+ return (eps_r2 * this->derivativeProbe(direction, p1, p2));
+}
+
+template
+KernelS SphericalDiffuse::exportKernelS_impl() const {
+ return std::bind(&SphericalDiffuse::kernelS,
+ *this,
+ std::placeholders::_1,
+ std::placeholders::_2);
+}
+
+template
+KernelD SphericalDiffuse::exportKernelD_impl() const {
+ return std::bind(&SphericalDiffuse::kernelD,
+ *this,
+ std::placeholders::_1,
+ std::placeholders::_2,
+ std::placeholders::_3);
+}
+
+template
+DerivativeProbe SphericalDiffuse::exportDerivativeProbe_impl()
+ const {
+ return std::bind(&SphericalDiffuse::derivativeProbe,
+ *this,
+ std::placeholders::_1,
+ std::placeholders::_2,
+ std::placeholders::_3);
+}
+
+template
+double SphericalDiffuse::singleLayer_impl(const Element & e,
+ double factor) const {
+ // Diagonal of S inside the cavity
+ double Sii_I = detail::diagonalSi(e.area(), factor);
+ // "Diagonal" of Coulomb singularity separation coefficient
+ double coulomb_coeff = this->coefficientCoulomb(e.center(), e.center());
+ // "Diagonal" of the image Green's function
+ double image = this->imagePotential(e.center(), e.center());
+ return (Sii_I / coulomb_coeff + image);
+}
+
+template
+double SphericalDiffuse::doubleLayer_impl(const Element & e,
+ double factor) const {
+ // Diagonal of S inside the cavity
+ double Sii_I = detail::diagonalSi(e.area(), factor);
+ // Diagonal of D inside the cavity
+ double Dii_I = detail::diagonalDi(e.area(), e.sphere().radius, factor);
+ // "Diagonal" of Coulomb singularity separation coefficient
+ double coulomb_coeff = this->coefficientCoulomb(e.center(), e.center());
+ // "Diagonal" of the directional derivative of the Coulomb singularity
+ // separation coefficient
+ double coeff_grad =
+ this->coefficientCoulombDerivative(e.normal(), e.center(), e.center()) /
+ std::pow(coulomb_coeff, 2);
+ // "Diagonal" of the directional derivative of the image Green's function
+ double image_grad =
+ this->imagePotentialDerivative(e.normal(), e.center(), e.center());
+
+ double eps_r2 = 0.0;
+ std::tie(eps_r2, std::ignore) = this->epsilon(e.center());
+
+ return (eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad));
+}
+
+template
+std::ostream & SphericalDiffuse::printObject(std::ostream & os) {
+ Eigen::IOFormat CleanFmt(Eigen::StreamPrecision, 0, ", ", "\n", "(", ")");
+ os << "Green's function type: spherical diffuse" << std::endl;
+ os << this->profile_ << std::endl;
+ os << "Sphere center = " << this->origin_.transpose().format(CleanFmt)
+ << std::endl;
+ os << "Angular momentum (Green's function) = " << this->maxLGreen_ << std::endl;
+ os << "Angular momentum (Coulomb coefficient) = " << this->maxLC_;
+ return os;
+}
+
+template
+void SphericalDiffuse::initSphericalDiffuse() {
+ using namespace detail;
+
+ // Parameters for the numerical solution of the radial differential equation
+ /*! 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();
+ /*! Time step between observer calls */
+ double step_ = 1.0e-2 * relative_width;
+
+ ProfileEvaluator eval_ =
+ std::bind(&ProfilePolicy::operator(), this->profile_, std::placeholders::_1);
+
+ zetaC_ =
+ RadialFunction(maxLC_, y_0_, y_infinity_, step_, eval_);
+ omegaC_ =
+ RadialFunction(maxLC_, y_infinity_, y_0_, step_, eval_);
+ zeta_.reserve(maxLGreen_ + 1);
+ omega_.reserve(maxLGreen_ + 1);
+ for (int L = 0; L <= maxLGreen_; ++L) {
+ zeta_.push_back(
+ RadialFunction(L, y_0_, y_infinity_, step_, eval_));
+ omega_.push_back(
+ RadialFunction(L, y_infinity_, y_0_, step_, eval_));
+ }
+}
+
+template
+double SphericalDiffuse::imagePotentialComponent_impl(
+ int L,
+ const Eigen::Vector3d & sp,
+ const Eigen::Vector3d & pp,
+ double Cr12) const {
+ Eigen::Vector3d sp_shift = sp + this->origin_;
+ Eigen::Vector3d pp_shift = pp + this->origin_;
+ double r1 = sp_shift.norm();
+ double r2 = pp_shift.norm();
+ double cos_gamma = sp_shift.dot(pp_shift) / (r1 * r2);
+ // Evaluate Legendre polynomial of order L
+ // First of all clean-up cos_gamma, Legendre polynomials
+ // are only defined for -1 <= x <= 1
+ if (utils::numericalZero(cos_gamma - 1))
+ cos_gamma = 1.0;
+ if (utils::numericalZero(cos_gamma + 1))
+ cos_gamma = -1.0;
+ double pl_x = Legendre::Pn(L, cos_gamma);
+
+ /* Zeta and Omega are now on a logarithmic scale We need to pass the
+ arguments in the correct scale and then use the chain rule to get
+ the proper derivative */
+ double y1 = std::log(r1);
+ double y2 = std::log(r2);
+
+ /* Sample zeta_[L] */
+ double zeta1 = 0.0, zeta2 = 0.0, d_zeta2 = 0.0;
+ /* Value of zeta_[L] at point with index 1 */
+ std::tie(zeta1, std::ignore) = zeta_[L](y1);
+ /* Value of zeta_[L] and its first derivative at point with index 2 */
+ std::tie(zeta2, d_zeta2) = zeta_[L](y2);
+
+ /* Sample omega_[L] */
+ double omega1 = 0.0, omega2 = 0.0, d_omega2 = 0.0;
+ /* Value of omega_[L] at point with index 1 */
+ std::tie(omega1, std::ignore) = omega_[L](y1);
+ /* Value of omega_[L] and its first derivative at point with index 2 */
+ std::tie(omega2, d_omega2) = omega_[L](y2);
+
+ double eps_r2 = 0.0;
+ std::tie(eps_r2, std::ignore) = this->profile_(pp_shift.norm());
+
+ /* Evaluation of the Wronskian and the denominator */
+ double denominator = (d_zeta2 - d_omega2) * r2 * eps_r2;
+
+ double gr12 = 0.0;
+ if (r1 < r2) {
+ gr12 = std::exp(zeta1 - zeta2) * (2 * L + 1) / denominator;
+ double f_L = r1 / r2;
+ for (int i = 1; i < L; ++i) {
+ f_L *= r1 / r2;
+ }
+ gr12 = (gr12 - f_L / (r2 * Cr12)) * pl_x;
+ } else {
+ gr12 = std::exp(omega1 - omega2) * (2 * L + 1) / denominator;
+ double f_L = r2 / r1;
+ for (int i = 1; i < L; ++i) {
+ f_L *= r2 / r1;
+ }
+ gr12 = (gr12 - f_L / (r1 * Cr12)) * pl_x;
+ }
+ return gr12;
+}
+
+template
+double SphericalDiffuse::coefficient_impl(
+ const Eigen::Vector3d & sp,
+ const Eigen::Vector3d & pp) const {
+ double r1 = (sp + this->origin_).norm();
+ double r2 = (pp + this->origin_).norm();
+
+ double y1 = std::log(r1);
+ double y2 = std::log(r2);
+
+ /* Sample zetaC_ */
+ double zeta1 = 0.0, zeta2 = 0.0, d_zeta2 = 0.0;
+ /* Value of zetaC_ at point with index 1 */
+ std::tie(zeta1, std::ignore) = zetaC_(y1);
+ /* Value of zetaC_ and its first derivative at point with index 2 */
+ std::tie(zeta2, d_zeta2) = zetaC_(y2);
+
+ /* Sample omegaC_ */
+ double omega1 = 0.0, omega2 = 0.0, d_omega2 = 0.0;
+ /* Value of omegaC_ at point with index 1 */
+ std::tie(omega1, std::ignore) = omegaC_(y1);
+ /* Value of omegaC_ and its first derivative at point with index 2 */
+ std::tie(omega2, d_omega2) = omegaC_(y2);
+
+ double tmp = 0.0, coeff = 0.0;
+ double eps_r2 = 0.0;
+ std::tie(eps_r2, std::ignore) = this->profile_(r2);
+
+ /* Evaluation of the Wronskian and the denominator */
+ double denominator = (d_zeta2 - d_omega2) * r2 * eps_r2;
+
+ if (r1 < r2) {
+ double f_L = r1 / r2;
+ for (int i = 1; i < maxLC_; ++i) {
+ f_L *= r1 / r2;
+ }
+ tmp = std::exp(zeta1 - zeta2) * (2 * maxLC_ + 1) / denominator;
+ coeff = f_L / (tmp * r2);
+ } else {
+ double f_L = r2 / r1;
+ for (int i = 1; i < maxLC_; ++i) {
+ f_L *= r2 / r1;
+ }
+ tmp = std::exp(omega1 - omega2) * (2 * maxLC_ + 1) / denominator;
+ coeff = f_L / (tmp * r1);
+ }
+
+ return coeff;
+}
+
+using dielectric_profile::OneLayerTanh;
+template class SphericalDiffuse;
+
+using dielectric_profile::OneLayerErf;
+template class SphericalDiffuse;
+
+using dielectric_profile::OneLayerLog;
+template class SphericalDiffuse;
+
+} // namespace green
+} // namespace pcm
diff --git a/src/green/SphericalDiffuse.hpp b/src/green/SphericalDiffuse.hpp
new file mode 100644
index 0000000..0567328
--- /dev/null
+++ b/src/green/SphericalDiffuse.hpp
@@ -0,0 +1,272 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 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
+
+/*! \file SphericalDiffuse.hpp */
+
+namespace pcm {
+namespace cavity {
+class Element;
+} // namespace cavity
+} // namespace pcm
+
+#include "GreenData.hpp"
+#include "GreensFunction.hpp"
+#include "InterfacesImpl.hpp"
+#include "dielectric_profile/OneLayerErf.hpp"
+#include "dielectric_profile/OneLayerLog.hpp"
+#include "dielectric_profile/OneLayerTanh.hpp"
+
+namespace pcm {
+namespace green {
+/*! \class SphericalDiffuse
+ * \brief Green's function for a diffuse interface with spherical symmetry
+ * \author Hui Cao, Ville Weijo, Luca Frediani and Roberto Di Remigio
+ * \date 2010-2015
+ * \tparam ProfilePolicy functional form of the diffuse layer
+ *
+ * The origin of the dielectric sphere can be changed by means of the constructor.
+ * The solution of the differential equation defining the Green's function is
+ * **always** performed assuming that the dielectric sphere is centered in the
+ * origin of the
+ * coordinate system. Whenever the public methods are invoked to "sample" the
+ * Green's function
+ * at a pair of points, a translation of the sampling points is performed first.
+ */
+template
+class SphericalDiffuse final : public GreensFunction {
+public:
+ /*! Constructor for a one-layer interface
+ * \param[in] e1 left-side dielectric constant
+ * \param[in] e2 right-side dielectric constant
+ * \param[in] w width of the interface layer
+ * \param[in] c center of the diffuse layer
+ * \param[in] o center of the sphere
+ * \param[in] l maximum value of angular momentum
+ */
+ SphericalDiffuse(double e1,
+ double e2,
+ double w,
+ double c,
+ const Eigen::Vector3d & o,
+ int l)
+ : GreensFunction(ProfilePolicy(e1, e2, w, c)),
+ origin_(o),
+ maxLGreen_(l),
+ maxLC_(2 * l) {
+ initSphericalDiffuse();
+ }
+
+ virtual double permittivity() const override final {
+ PCMSOLVER_ERROR("permittivity() only implemented for uniform dielectrics");
+ }
+
+ /*! \brief Returns Coulomb singularity separation coefficient
+ * \param[in] source location of the source charge
+ * \param[in] probe location of the probe charge
+ */
+ double coefficientCoulomb(const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const;
+
+ /*! \brief Returns singular part of the Green's function
+ * \param[in] source location of the source charge
+ * \param[in] probe location of the probe charge
+ */
+ double Coulomb(const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const;
+
+ /*! \brief Returns non-singular part of the Green's function (image potential)
+ * \param[in] source location of the source charge
+ * \param[in] probe location of the probe charge
+ */
+ double imagePotential(const Eigen::Vector3d & source,
+ const Eigen::Vector3d & probe) const;
+
+ /*! Returns value of the directional derivative of the
+ * Coulomb singularity separation coefficient for the pair of points p1, p2:
+ * \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+ *\mathbf{n}_{\mathbf{p}_2}\f$
+ * Notice that this method returns the directional derivative with respect
+ * to the probe point, thus assuming that the direction is relative to that
+ *point.
+ *
+ * \param[in] direction the direction
+ * \param[in] p1 first point
+ * \param[in] p2 second point
+ */
+ double coefficientCoulombDerivative(const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const;
+
+ /*! Returns value of the directional derivative of the
+ * singular part of the Greens's function for the pair of points p1, p2:
+ * \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+ *\mathbf{n}_{\mathbf{p}_2}\f$
+ * Notice that this method returns the directional derivative with respect
+ * to the probe point, thus assuming that the direction is relative to that
+ *point.
+ *
+ * \param[in] direction the direction
+ * \param[in] p1 first point
+ * \param[in] p2 second point
+ */
+ double CoulombDerivative(const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const;
+
+ /*! Returns value of the directional derivative of the
+ * non-singular part (image potential) of the Greens's function for the pair of
+ *points p1, p2:
+ * \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+ *\mathbf{n}_{\mathbf{p}_2}\f$
+ * Notice that this method returns the directional derivative with respect
+ * to the probe point, thus assuming that the direction is relative to that
+ *point.
+ *
+ * \param[in] direction the direction
+ * \param[in] p1 first point
+ * \param[in] p2 second point
+ */
+ double imagePotentialDerivative(const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const;
+
+ /*! Handle to the dielectric profile evaluation */
+ std::tuple epsilon(const Eigen::Vector3d & point) const;
+ void toFile(const std::string & prefix = "");
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+private:
+ /*! Evaluates the Green's function given a pair of points
+ * \param[in] sp the source point
+ * \param[in] pp the probe point
+ *
+ * \note This takes care of the origin shift
+ */
+ virtual Stencil operator()(Stencil * sp, Stencil * pp) const override;
+
+ /*! Returns value of the kernel of the \f$\mathcal{D}\f$ integral operator for the
+ * pair of points p1, p2:
+ * \f$ [\boldsymbol{\varepsilon}\nabla_{\mathbf{p_2}}G(\mathbf{p}_1,
+ * \mathbf{p}_2)]\cdot \mathbf{n}_{\mathbf{p}_2}\f$
+ * To obtain the kernel of the \f$\mathcal{D}^\dagger\f$ operator call this
+ * methods with \f$\mathbf{p}_1\f$
+ * and \f$\mathbf{p}_2\f$ exchanged and with \f$\mathbf{n}_{\mathbf{p}_2} =
+ * \mathbf{n}_{\mathbf{p}_1}\f$
+ * \param[in] direction the direction
+ * \param[in] p1 first point
+ * \param[in] p2 second point
+ */
+ virtual double kernelD_impl(const Eigen::Vector3d & direction,
+ const Eigen::Vector3d & p1,
+ const Eigen::Vector3d & p2) const override;
+
+ virtual KernelS exportKernelS_impl() const override;
+
+ virtual KernelD exportKernelD_impl() const override;
+
+ virtual DerivativeProbe exportDerivativeProbe_impl() const override;
+
+ virtual double singleLayer_impl(const Element & e, double factor) const override;
+
+ virtual double doubleLayer_impl(const Element & e, double factor) const override;
+
+ virtual std::ostream & printObject(std::ostream & os) override;
+ /*! This calculates all the components needed to evaluate the Green's function */
+ void initSphericalDiffuse();
+
+ /*! Center of the dielectric sphere */
+ Eigen::Vector3d origin_;
+
+ /*! @{ Parameters and functions for the calculation of the Green's function,
+ * including Coulomb singularity */
+ /*! Maximum angular momentum in the final summation over Legendre polynomials to
+ * obtain G */
+ int maxLGreen_;
+
+ /*! \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_;
+
+ /*! \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_;
+
+ /*! \brief Returns L-th component of the radial part of the Green's function
+ * \param[in] L angular momentum
+ * \param[in] sp source point
+ * \param[in] pp probe point
+ * \param[in] Cr12 Coulomb singularity separation coefficient
+ * \note This function shifts the given source and probe points by the location
+ * of
+ * the
+ * dielectric sphere.
+ */
+ double imagePotentialComponent_impl(int L,
+ const Eigen::Vector3d & sp,
+ const Eigen::Vector3d & pp,
+ double Cr12) const;
+ /*! @}*/
+ /*! @{ Parameters and functions for the calculation of the Coulomb singularity
+ * separation coefficient */
+ /*! Maximum angular momentum to obtain C(r, r'), needed to separate the Coulomb
+ * singularity */
+ int maxLC_; // = 2 * maxLGreen_;
+
+ /*! \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_;
+
+ /*! \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_;
+
+ /*! \brief Returns coefficient for the separation of the Coulomb singularity
+ * \param[in] sp first point
+ * \param[in] pp second point
+ * \note This function shifts the given source and probe points by the location
+ * of
+ * the
+ * dielectric sphere.
+ */
+ double coefficient_impl(const Eigen::Vector3d & sp,
+ const Eigen::Vector3d & pp) const;
+ /*! @}*/
+};
+
+template
+IGreensFunction * createSphericalDiffuse(const GreenData & data) {
+ return new SphericalDiffuse(
+ data.epsilon1, data.epsilon2, data.width, data.center, data.origin, data.maxL);
+}
+} // namespace green
+} // namespace pcm
diff --git a/src/interface/Meddle.cpp b/src/interface/Meddle.cpp
index f8add79..ec00da5 100644
--- a/src/interface/Meddle.cpp
+++ b/src/interface/Meddle.cpp
@@ -442,7 +442,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 bc75099..f58e200 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -22,6 +22,7 @@ list(APPEND headers_list
Molecule.hpp
PhysicalConstants.hpp
QuadratureRules.hpp
+ RungeKutta.hpp
Solvent.hpp
Sphere.hpp
SplineFunction.hpp
diff --git a/src/utils/MathUtils.hpp b/src/utils/MathUtils.hpp
index 24c8042..03b0538 100644
--- a/src/utils/MathUtils.hpp
+++ b/src/utils/MathUtils.hpp
@@ -44,48 +44,6 @@
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
@@ -112,110 +70,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,
- int cavitySize,
- int 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 (int a = 0; a < cavitySize; ++a) {
- for (int 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
*
@@ -245,11 +100,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
*/
@@ -264,69 +119,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 bf21fa9..dab578f 100644
--- a/src/utils/Molecule.cpp
+++ b/src/utils/Molecule.cpp
@@ -56,7 +56,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,
@@ -74,7 +74,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,
@@ -107,7 +107,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 8c9a9ec..6a3628e 100644
--- a/src/utils/Molecule.hpp
+++ b/src/utils/Molecule.hpp
@@ -46,6 +46,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",
@@ -90,7 +91,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 0000000..4f0fb36
--- /dev/null
+++ b/src/utils/RungeKutta.hpp
@@ -0,0 +1,195 @@
+/*
+ * 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
+
+/*! \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 c35a405..0000000
--- 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 8759168..0000000
--- 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 79284c3..c9eba69 100644
--- a/src/utils/SplineFunction.hpp
+++ b/src/utils/SplineFunction.hpp
@@ -24,19 +24,22 @@
#pragma once
#include
+#include
#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;
@@ -84,3 +87,43 @@ class SplineFunction final {
.transpose();
}
};
+
+/*! \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 efb6098..090b538 100644
--- a/src/utils/Symmetry.cpp
+++ b/src/utils/Symmetry.cpp
@@ -23,25 +23,91 @@
#include "Symmetry.hpp"
+#include
+
#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,
+ int cavitySize,
+ int 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 (int a = 0; a < cavitySize; ++a) {
+ for (int 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 4d6aa15..e26bc26 100644
--- a/src/utils/Symmetry.hpp
+++ b/src/utils/Symmetry.hpp
@@ -24,57 +24,134 @@
#pragma once
#include
+#include
#include
+#include
-/*! \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 final {
+public:
+ 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}}) {}
+ 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];
+ }
+
+ int nrGenerators() const { return nrGenerators_; }
+ int nrIrrep() const { return nrIrrep_; }
+ int generators(size_t i) const { return generators_[i]; }
-class Symmetry {
private:
- /*!
- * Number of generators
- */
+ /*! Number of generators */
int nrGenerators_;
- /*!
- * Generators
- */
- int generators_[3];
- /*!
- * Number of irreps
- */
+ /*! Number of irreps */
int nrIrrep_;
+ /*! Generators */
+ std::array generators_;
+};
-public:
- /*! \brief Default constructor sets up C1 point group
- */
- Symmetry() : nrGenerators_(0) {
- 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(const Symmetry & other)
- : nrGenerators_(other.nrGenerators_), nrIrrep_(other.nrIrrep_) {
- std::copy(other.generators_, other.generators_ + nrGenerators_, generators_);
+/*! \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,
+ int cavitySize,
+ int 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];
}
- ~Symmetry() {}
- int nrGenerators() const { return nrGenerators_; }
- int generators(int i) const { return generators_[i]; }
- int nrIrrep() const { return nrIrrep_; }
-};
+ return parity;
+}
-/*! Builds Symmetry object.
+/*! \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 38e4b3d..bb1e8e4 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);
};