Get some work done on the Rob-green_spherical branch

Author: Roberto Di Remigio

src/bi_operators/CollocationIntegrator.hpp

@@ -2,6 +2,7 @@
 #define COLLOCATIONINTEGRATOR_HPP
 
 #include <iosfwd>
+#include <tuple>
 
 #include "Config.hpp"
 
@@ -54,11 +55,10 @@ struct CollocationIntegrator
 	    double area = e.area();
 	    return (factor_ * std::sqrt(4 * M_PI / area) * epsInv);
     }
-    double computeD(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const Element & e) const {
-        double epsInv = 1.0 / gf.epsilon();
+    double computeD(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & e) const {
 	    double area = e.area();
 	    double radius = e.sphere().radius();
-        return (-factor_ * std::sqrt(M_PI/ area) * (1.0 / radius) * epsInv);
+        return (-factor_ * std::sqrt(M_PI/ area) * (1.0 / radius));
     }
 
     double computeS(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & /* e */) const {
@@ -102,7 +102,10 @@ struct CollocationIntegrator
         // "Diagonal" of the directional derivative of the image Green's function
         double image_grad = gf.imagePotentialDerivative(e.normal(), e.center(), e.center());
 
-        return (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
+        double eps_r2 = 0.0;
+        std::tie(eps_r2, std::ignore) = gf.epsilon(e.center());
+
+        return eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
     }
 
     /// Scaling factor for the collocation formulas

src/green/SphericalDiffuse.hpp

@@ -157,16 +157,13 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
     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_;
-        double r1  = source_shifted.norm();
-        double r2  = probe_shifted.norm();
-        double cos_gamma = source_shifted.dot(probe_shifted) / (r1 * r2);
 
         // Obtain coefficient for the separation of the Coulomb singularity
-        double Cr12 = this->coefficient(r1, r2);
+        double Cr12 = this->coefficient(source_shifted.norm(), probe_shifted.norm());
 
         double gr12 = 0.0;
         for (int L = 0; L <= maxLGreen_; ++L) {
-            gr12 += this->functionSummation(L, r1, r2, cos_gamma, Cr12);
+            gr12 += this->functionSummation(L, source_shifted, probe_shifted, Cr12);
         }
 
         return gr12;
@@ -217,7 +214,9 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
                 p2, p1, direction, this->delta_);
     }
     /*! Handle to the dielectric profile evaluation */
-    void epsilon(double & v, double & d, double point) const { std::tie(v, d) = this->profile_(point); }
+    std::tuple<double, double> epsilon(const Eigen::Vector3d & point) const {
+        return this->profile_((point - this->origin_).norm());
+    }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW /* See http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html */
 private:
     using RadialFunction = interfaces::RadialFunction;
@@ -233,18 +232,15 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
         Eigen::Map<Eigen::Matrix<double, 3, 1> > p1(sp), p2(pp);
         Eigen::Vector3d source = p1 - this->origin_;
         Eigen::Vector3d probe  = p2 - this->origin_;
-        double r1  = source.norm();
-        double r2  = probe.norm();
-        double r12 = (source - probe).norm();
-        double cos_gamma = source.dot(probe) / (r1 * r2);
 
         // Obtain coefficient for the separation of the Coulomb singularity
-        double Cr12 = this->coefficient(r1, r2);
+        double Cr12 = this->coefficient(source.norm(), probe.norm());
 
         double gr12 = 0.0;
         for (int L = 0; L <= maxLGreen_; ++L) {
-            gr12 += this->functionSummation(L, r1, r2, cos_gamma, Cr12);
+            gr12 += this->functionSummation(L, source, probe, Cr12);
         }
+        double r12 = (source - probe).norm();
 
         return (1.0 / (Cr12 * r12) + gr12);
     }
@@ -335,13 +331,16 @@ 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] r1 first point
-     *  \param[in] r2 second point
-     *  \param[in] cos_gamma cosine of the angle in between first and second point
+     *  \param[in] sp source point, shifted by origin_
+     *  \param[in] pp probe point, shifted by origin_
      *  \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.
      */
-    double functionSummation(int L, double r1, double r2, double cos_gamma, double Cr12) const {
-        double gr12 = 0.0;
+    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);
         // Evaluate Legendre polynomial of order L
         // First of all clean-up cos_gamma, Legendre polynomials
         // are only defined for -1 <= x <= 1
@@ -369,6 +368,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
 
         double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
 
+        double gr12 = 0.0;
         if (r1 < r2) {
             gr12 = std::exp(zeta1 - zeta2) * (2*L +1) / denominator;
             gr12 = (gr12 - std::pow(r1/r2, L) / (r2 * Cr12) ) * pl_x ;

src/solver/PCMSolver.hpp

@@ -54,10 +54,10 @@ class PCMSolver
     PCMSolver(IGreensFunction * gfInside, IGreensFunction * gfOutside)
             : greenInside_(gfInside), greenOutside_(gfOutside), allocated_(true) {}
     virtual ~PCMSolver() {
-        if (allocated_) {
-            delete greenInside_;
-            delete greenOutside_;
-        }
+      //if (allocated_) {
+      //    delete greenInside_;
+      //    delete greenOutside_;
+      //}
     }
 
     IGreensFunction * greenInside() const {

tests/bi_operators/CMakeLists.txt

@@ -1,6 +1,10 @@
 set(list_of_libraries green cavity pedra utils)
 
+
 list(APPEND external_libraries ${ZLIB_LIBRARIES} ${Boost_SYSTEM_LIBRARY} ${Boost_TIMER_LIBRARY} ${Boost_CHRONO_LIBRARY})
+add_executable(update_reference_files.x update_reference_files.cpp)
+target_link_libraries(update_reference_files.x "${list_of_libraries}" "${external_libraries}")
+
 add_boosttest(bi_operators_collocation "${list_of_libraries}" "${external_libraries}")
 add_reference(vacuum_S_collocation.npy ${CMAKE_CURRENT_BINARY_DIR})
 add_reference(vacuum_D_collocation.npy ${CMAKE_CURRENT_BINARY_DIR})

tests/bi_operators/bi_operators_collocation.cpp

@@ -81,12 +81,6 @@ BOOST_FIXTURE_TEST_CASE(vacuum, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	S_results(i) = gf.diagonalS(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    unsigned int dim = static_cast<unsigned int>(cavity.size());
-    const unsigned int shape[] = {dim};
-    cnpy::npy_save("vacuum_S_collocation.npy", S_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_S_ref = cnpy::npy_load("vacuum_S_collocation.npy");
     int dim_read = raw_S_ref.shape[0];
     Eigen::VectorXd S_reference = Eigen::VectorXd::Zero(dim_read);
@@ -104,10 +98,6 @@ BOOST_FIXTURE_TEST_CASE(vacuum, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	D_results(i) = gf.diagonalD(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    cnpy::npy_save("vacuum_D_collocation.npy", D_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_D_ref = cnpy::npy_load("vacuum_D_collocation.npy");
     dim_read = raw_D_ref.shape[0];
     Eigen::VectorXd D_reference = Eigen::VectorXd::Zero(dim_read);
@@ -134,12 +124,6 @@ BOOST_FIXTURE_TEST_CASE(uniformdielectric, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	S_results(i) = gf.diagonalS(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    unsigned int dim = static_cast<unsigned int>(cavity.size());
-    const unsigned int shape[] = {dim};
-    cnpy::npy_save("uniformdielectric_S_collocation.npy", S_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_S_ref = cnpy::npy_load("uniformdielectric_S_collocation.npy");
     int dim_read = raw_S_ref.shape[0];
     Eigen::VectorXd S_reference = Eigen::VectorXd::Zero(dim_read);
@@ -157,10 +141,6 @@ BOOST_FIXTURE_TEST_CASE(uniformdielectric, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	D_results(i) = gf.diagonalD(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    cnpy::npy_save("uniformdielectric_D_collocation.npy", D_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_D_ref = cnpy::npy_load("uniformdielectric_D_collocation.npy");
     dim_read = raw_D_ref.shape[0];
     Eigen::VectorXd D_reference = Eigen::VectorXd::Zero(dim_read);
@@ -186,12 +166,6 @@ BOOST_FIXTURE_TEST_CASE(tanhsphericaldiffuse, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	    S_results(i) = gf.diagonalS(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    unsigned int dim = static_cast<unsigned int>(cavity.size());
-    const unsigned int shape[] = {dim};
-    cnpy::npy_save("tanhsphericaldiffuse_S_collocation.npy", S_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_S_ref = cnpy::npy_load("tanhsphericaldiffuse_S_collocation.npy");
     int dim_read = raw_S_ref.shape[0];
     Eigen::VectorXd S_reference = Eigen::VectorXd::Zero(dim_read);
@@ -209,10 +183,6 @@ BOOST_FIXTURE_TEST_CASE(tanhsphericaldiffuse, CollocationIntegratorTest)
     for (int i = 0; i < cavity.size(); ++i) {
 	    D_results(i) = gf.diagonalD(cavity.elements(i));
     }
-    // In case you need to update the reference files...
-    /*
-    cnpy::npy_save("tanhsphericaldiffuse_D_collocation.npy", D_results.data(), shape, 1, "w", false);
-    */
     cnpy::NpyArray raw_D_ref = cnpy::npy_load("tanhsphericaldiffuse_D_collocation.npy");
     dim_read = raw_D_ref.shape[0];
     Eigen::VectorXd D_reference = Eigen::VectorXd::Zero(dim_read);

tests/bi_operators/tanhsphericaldiffuse_D_collocation.npy

@@ -0,0 +1,131 @@
+/* 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 */
+
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <limits>
+
+#include "Config.hpp"
+
+#include <Eigen/Dense>
+
+#include "cnpyPimpl.hpp"
+#include "CollocationIntegrator.hpp"
+#include "DerivativeTypes.hpp"
+#include "Element.hpp"
+#include "GePolCavity.hpp"
+#include "PhysicalConstants.hpp"
+#include "SphericalDiffuse.hpp"
+#include "TestingMolecules.hpp"
+#include "UniformDielectric.hpp"
+#include "Vacuum.hpp"
+
+void save_vacuum();
+void save_uniform_dielectric();
+void save_tanh_spherical_diffuse();
+
+int main()
+{
+    save_vacuum();
+    save_uniform_dielectric();
+    save_tanh_spherical_diffuse();
+
+    return EXIT_SUCCESS;
+}
+
+void save_vacuum() {
+    Molecule molec = dummy<0>(1.44 / convertBohrToAngstrom);
+    double area = 10.0;
+    GePolCavity cavity(molec, area, 0.0, 100.0);
+
+    unsigned int dim = static_cast<unsigned int>(cavity.size());
+    const unsigned int shape[] = {dim};
+
+    Vacuum<AD_directional, CollocationIntegrator<AD_directional, Uniform> > gf;
+
+    Eigen::VectorXd S_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        S_results(i) = gf.diagonalS(cavity.elements(i));
+    }
+    cnpy::npy_save("vacuum_S_collocation.npy", S_results.data(), shape, 1, "w", false);
+
+    Eigen::VectorXd D_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        D_results(i) = gf.diagonalD(cavity.elements(i));
+    }
+    cnpy::npy_save("vacuum_D_collocation.npy", D_results.data(), shape, 1, "w", false);
+}
+
+void save_uniform_dielectric() {
+    double epsilon = 80.0;
+    Molecule molec = dummy<0>(1.44 / convertBohrToAngstrom);
+    double area = 10.0;
+    GePolCavity cavity(molec, area, 0.0, 100.0);
+
+    unsigned int dim = static_cast<unsigned int>(cavity.size());
+    const unsigned int shape[] = {dim};
+
+    UniformDielectric<AD_directional, CollocationIntegrator<AD_directional, Uniform> > gf(epsilon);
+
+    Eigen::VectorXd S_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        S_results(i) = gf.diagonalS(cavity.elements(i));
+    }
+    cnpy::npy_save("uniformdielectric_S_collocation.npy", S_results.data(), shape, 1, "w", false);
+
+    Eigen::VectorXd D_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        D_results(i) = gf.diagonalD(cavity.elements(i));
+    }
+    cnpy::npy_save("uniformdielectric_D_collocation.npy", D_results.data(), shape, 1, "w", false);
+}
+
+void save_tanh_spherical_diffuse() {
+    double epsilon = 80.0;
+    double width = 5.0;
+    double sphereRadius = 100.0;
+    Molecule molec = dummy<0>(1.44 / convertBohrToAngstrom);
+    double area = 10.0;
+    GePolCavity cavity(molec, area, 0.0, 100.0);
+
+    unsigned int dim = static_cast<unsigned int>(cavity.size());
+    const unsigned int shape[] = {dim};
+
+    SphericalDiffuse<CollocationIntegrator<Numerical, TanhDiffuse>, TanhDiffuse> gf(epsilon, epsilon, width, sphereRadius, Eigen::Vector3d::Zero());
+
+    Eigen::VectorXd S_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        S_results(i) = gf.diagonalS(cavity.elements(i));
+    }
+    cnpy::npy_save("tanhsphericaldiffuse_S_collocation.npy", S_results.data(), shape, 1, "w", false);
+
+    Eigen::VectorXd D_results = Eigen::VectorXd::Zero(cavity.size());
+    for (int i = 0; i < cavity.size(); ++i) {
+        D_results(i) = gf.diagonalD(cavity.elements(i));
+    }
+    cnpy::npy_save("tanhsphericaldiffuse_D_collocation.npy", D_results.data(), shape, 1, "w", false);
+}