Merge branch 'Rob-green_spherical' of gitlab.com:PCMSolver/pcmsolver into Rob-green_spherical

Author: Roberto Di Remigio

src/bin/plot_green_spherical.cpp

@@ -14,7 +14,6 @@
 #include "UniformDielectric.hpp"
 
 std::vector<double> generate_grid(double, double, size_t);
-std::string header();
 void case1();
 void case2();
 void case3();
@@ -40,59 +39,48 @@ std::vector<double> generate_grid(double min_val, double max_val, size_t nPoints
     return grid;
 }
 
-std::string header()
-{
-    std::stringstream tmp;
-    tmp << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
-    tmp << "gf_value" << '\t' << "image" << '\t' << "Coulomb" << '\t';
-    tmp << "derivativeProbe" << '\t' << "der_image" << '\t' << "der_Coulomb" << std::endl;
-    return tmp.str();
-}
-
 void case1()
 {
-    size_t nPoints = 100;
+    size_t nPoints = 200;
     double epsInside = 80.0;
     double epsOutside = 2.0;
     Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
     double sphereRadius = 100.0;
-    double width = 10.0;
+    double width = 9.0;
 
-    Eigen::Vector3d source;
-    source << 4.0, 0.0, 85.0;
+    Eigen::Vector3d source = Eigen::Vector3d::Zero();
     Eigen::Vector3d probe = Eigen::Vector3d::Zero();
 
-    double x_min = 10.0;
-    double x_max = 300.0;
-    std::vector<double> x = generate_grid(x_min, x_max, nPoints);
-    double z_min = 10.0;
-    double z_max = 300.0;
+    double z_min = -20.0;
+    double z_max = 20.0;
     std::vector<double> z = generate_grid(z_min, z_max, nPoints);
+    double delta = 0.1;
 
     std::ofstream out;
 
     SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
     LOG(gf);
 
     out.open("gf_spherical_CASE1.dat");
-    out << header();
     out.precision(16);
     for (size_t i = 0; i < nPoints; ++i) {
+	source << 0.0, 0.0, z[i];
         for (size_t j = 0; j < nPoints; ++j) {
-            probe << x[i], 0.0, z[j];
-            out << '\t' << x[i]
+            probe << 0.0, 0.0, (z[j] + delta);
+            out << '\t' << i
+		<< '\t' << j
+		<< '\t' << z[i]
                 << '\t' << z[j]
-                << '\t' << gf.kernelS(source, probe)
-                << '\t' << gf.imagePotential(source, probe)
                 << '\t' << gf.Coulomb(source, probe)
-                << '\t' << gf.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe)
-                << '\t' << gf.CoulombDerivative(Eigen::Vector3d::UnitZ(), source, probe)
-                << '\t' << gf.imagePotentialDerivative(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
+                << '\t' << gf.imagePotential(source, probe)
+                << '\t' << gf.kernelS(source, probe)
+		<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe)) << std::endl;
         }
     }
     out.close();
     LOG_TIME;
 
+    /*
     UniformDielectric<AD_directional, CollocationIntegrator> gf_inside(epsInside);
     out.open("gf_uniform_inside_CASE1.dat");
     out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
@@ -124,4 +112,5 @@ void case1()
         }
     }
     out.close();
+    */
 }

src/green/SphericalDiffuse.hpp

@@ -101,8 +101,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
                               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->profile_((p2 - this->origin_).norm());
+        std::tie(eps_r2, std::ignore) = this->epsilon(p2);
 
         return (eps_r2 * this->derivativeProbe(direction, p1, p2));
     }
@@ -130,40 +129,30 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
      *  \param[in] probe location of the probe charge
      */
     double coefficientCoulomb(const Eigen::Vector3d & source, const Eigen::Vector3d & probe) const {
-        double r1  = (source - this->origin_).norm();
-        double r2  = (probe  - this->origin_).norm();
-
         // Obtain coefficient for the separation of the Coulomb singularity
-        return this->coefficient(r1, r2);
+        return this->coefficient_impl(source, probe);
     }
     /*! \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 {
-        Eigen::Vector3d source_shifted = source - this->origin_;
-        Eigen::Vector3d probe_shifted  = probe  - this->origin_;
-        double r1  = source_shifted.norm();
-        double r2  = probe_shifted.norm();
-        double r12 = (source_shifted - probe_shifted).norm();
+        double r12 = (source - probe).norm();
 
         // Obtain coefficient for the separation of the Coulomb singularity
-        return (1.0 / (this->coefficient(r1, r2) * r12));
+        return (1.0 / (this->coefficient_impl(source, probe) * r12));
     }
     /*! \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 {
-        Eigen::Vector3d source_shifted = source - this->origin_;
-        Eigen::Vector3d probe_shifted  = probe  - this->origin_;
-
         // Obtain coefficient for the separation of the Coulomb singularity
-        double Cr12 = this->coefficient(source_shifted.norm(), probe_shifted.norm());
+        double Cr12 = this->coefficient_impl(source, probe);
 
         double gr12 = 0.0;
-        for (int L = 0; L <= maxLGreen_; ++L) {
-            gr12 += this->functionSummation(L, source_shifted, probe_shifted, Cr12);
+        for (int L = 1; L <= maxLGreen_; ++L) {
+            gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
         }
 
         return gr12;
@@ -215,7 +204,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
     }
     /*! Handle to the dielectric profile evaluation */
     std::tuple<double, double> epsilon(const Eigen::Vector3d & point) const {
-        return this->profile_((point - this->origin_).norm());
+        return this->profile_((point + this->origin_).norm());
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW /* See http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html */
 private:
@@ -229,16 +218,14 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
     virtual Numerical operator()(Numerical * sp, Numerical * pp) const
     {
         // Transfer raw arrays to Eigen vectors using the Map type
-        Eigen::Map<Eigen::Matrix<double, 3, 1> > p1(sp), p2(pp);
-        Eigen::Vector3d source = p1 - this->origin_;
-        Eigen::Vector3d probe  = p2 - this->origin_;
+        Eigen::Map<Eigen::Matrix<double, 3, 1> > source(sp), probe(pp);
 
         // Obtain coefficient for the separation of the Coulomb singularity
-        double Cr12 = this->coefficient(source.norm(), probe.norm());
+        double Cr12 = this->coefficient_impl(source, probe);
 
         double gr12 = 0.0;
-        for (int L = 0; L <= maxLGreen_; ++L) {
-            gr12 += this->functionSummation(L, source, probe, Cr12);
+        for (int L = 1; L <= maxLGreen_; ++L) {
+            gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
         }
         double r12 = (source - probe).norm();
 
@@ -333,16 +320,18 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
     std::vector<RadialFunction> 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, shifted by origin_
-     *  \param[in] pp probe point, shifted by origin_
+     *  \param[in] sp source point
+     *  \param[in] pp probe point
      *  \param[in] Cr12 Coulomb singularity separation coefficient
-     *  \note This function expects that the source and probe points have been correctly shifted
-     *  to account for the position of the origin of the dielectric sphere.
+     *  \note This function shifts the given source and probe points by the location of the
+     *  dielectric sphere.
      */
-    double functionSummation(int L, const Eigen::Vector3d & sp, const Eigen::Vector3d & pp, double Cr12) const {
-        double r1 = sp.norm();
-        double r2 = pp.norm();
-        double cos_gamma = sp.dot(pp) / (r1 * r2);
+    double 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
@@ -366,7 +355,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
         double d_omega2 = linearInterpolation(r2, omega_[L][0], omega_[L][2]);
 
         double eps_r2 = 0.0;
-        std::tie(eps_r2, std::ignore) = this->profile_(r2);
+        std::tie(eps_r2, std::ignore) = this->profile_(pp_shift.norm());
 
         double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
 
@@ -395,10 +384,14 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
      */
     RadialFunction omegaC_;
     /*! \brief Returns coefficient for the separation of the Coulomb singularity
-     *  \param[in] r1 first point
-     *  \param[in] r2 second point
+     *  \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(double r1, double r2) const {
+    double coefficient_impl(const Eigen::Vector3d & sp, const Eigen::Vector3d & pp) const {
+	double r1 = (sp + this->origin_).norm();
+	double r2 = (pp + this->origin_).norm();
         /* Value of zetaC_ at point with index 1 */
         double zeta1  = linearInterpolation(r1, zetaC_[0], zetaC_[1]);
         /* Value of zetaC_ at point with index 2 */

src/green/dielectric_profile/OneLayerErf.hpp

@@ -65,7 +65,7 @@ class OneLayerErf
      *  \param[in] point where to evaluate the derivative
      */
     double derivative(double point) const {
-        double factor = (epsilon1_ - epsilon2_) / (width_ * std::sqrt(M_PI));
+        double factor = (epsilon2_ - epsilon1_) / (width_ * std::sqrt(M_PI));
         double t = (point - center_) / width_;
         double val = std::exp(-std::pow(t, 2));
         return (factor * val); //first derivative of epsilon(r)

src/green/dielectric_profile/OneLayerTanh.hpp

@@ -64,7 +64,7 @@ class OneLayerTanh
      *  \param[in] point where to evaluate the derivative
      */
     double derivative(double point) const {
-        double factor = (epsilon1_ - epsilon2_) / (2.0 * width_);
+        double factor = (epsilon2_ - epsilon1_) / (2.0 * width_);
         double tanh_r = std::tanh((point - center_) / width_);
         return (factor * (1 - std::pow(tanh_r, 2))); //first derivative of epsilon(r)
     }

src/utils/SplineFunction.hpp

@@ -0,0 +1,78 @@
+/* pcmsolver_copyright_start */
+/*
+ *     PCMSolver, an API for the Polarizable Continuum Model
+ *     Copyright (C) 2013 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 <http://www.gnu.org/licenses/>.
+ *
+ *     For information on the complete list of contributors to the
+ *     PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
+ */
+/* pcmsolver_copyright_end */
+
+#ifndef SPLINEFUNCTION_HPP
+#define SPLINEFUNCTION_HPP
+
+#include "Config.hpp"
+
+#include <Eigen/Core>
+#include <unsupported/Eigen/Splines>
+
+/*! \file SplineFunction.hpp
+ *  \class SplineFunction
+ *  \brief Spline interpolation of a function
+ *  \author Roberto Di Remigio
+ *  \date 2015
+ *
+ *  Taken from StackOverflow http://stackoverflow.com/a/29825204/2528668
+ */
+   
+/*! Scale abscissa value to [0, 1] */
+inline double scale(double x) {
+	return (x - xMin_) / (xMax_ - xMin_);
+}
+
+/*! Scale abscissa values to [0, 1] */
+inline Eigen::RowVectorXd scale(const Eigen::VectorXd & x_vec) {
+    return x_vec.unaryExpr([](double x) -> double { return scale(x); }).transpose();
+}
+
+class SplineFunction
+{
+private:
+    typedef Eigen::Spline<double, 1, 3> SplineFit;
+public:
+    /*! \param[in] x vector with abscissa values 
+     *  \param[in] y vector with function values
+     */
+    SplineFunction(const Eigen::VectorXd & x, const Eigen::VectorXd & y) 
+	    : xMin_(x.minCoeff()), xMax_(x.maxCoeff()),
+	    spline_(Eigen::SplineFitting<SplineFit>::Interpolate(y.transpose(), 
+		      std::min<int>(x.rows() - 1, 3), scale(x))
+	           )
+	{}
+    /*! \brief Evaluate spline at given point
+     *  \param[in] x evaluation point
+     */
+    double operator()(double x) const {
+	return spline_(scale(x))(0);
+    }
+private:
+    double xMin_;
+    double xMax_;
+    SplineFit spline_;
+};
+#endif // SPLINEFUNCTION_HPP