Users can now select the scaling factor for the diagonal of the collocation matrices

Author: Roberto Di Remigio

CHANGELOG.md

@@ -5,9 +5,12 @@
 - A runtime check to ensure that all atoms have a nonzero radius.
 API kills program execution if this is the case.
 - An API function to retrieve the areas/weights of the cavity finite elements.
-The values in the returned array are in Bohr^2. Addresses feature request from @shofener (Issue #13)
+The values in the returned array are in Bohr^2. Addresses a feature request from @shofener (Issue #13)
 
 ### Changed
+- Boundary integral operators classes learnt to accept a scaling factor for the diagonal
+elements of the approximate collocation matrices. The change is reflected in the
+Green's funtion classes and in the input parsing. Addresses a feature request from @shofener (Issue #16)
 - GePolCavity learnt to print also the list of spheres used to generate
 the cavity.
 - Different internal handling of conversion factors from Bohr to Angstrom.

doc/users/input.rst

@@ -241,10 +241,28 @@ Medium section keywords
      k}`
 
      * **Type**: double
-     * **Valide values**: :math:`k > 0.0`
+     * **Valid values**: :math:`k > 0.0`
      * **Valid for**: CPCM solver
      * **Default**: 0.0
 
+   DiagonalIntegrator
+     Type of integrator for the diagonal of the boundary integral operators
+
+     * **Type**: string
+     * **Valid values**: COLLOCATION
+     * **Valid for**: IEFPCM, CPCM
+     * **Default**: COLLOCATION
+     * **Notes**: in future releases we will add PURISIMA and NUMERICAL as options
+
+   DiagonalScaling
+     Scaling factor for diagonal of collocation matrices
+
+     * **Type**: double
+     * **Valid values**: :math:`f > 0.0`
+     * **Valid for**: IEFPCM, CPCM
+     * **Default**: 1.07
+     * **Notes**: values commonly used in the literature are 1.07 and 1.0694
+
    ProbeRadius
      Radius of the spherical probe approximating a solvent molecule. Used for
      generating the solvent-excluded surface (SES) or an approximation of it.

src/bi_operators/CollocationIntegrator.hpp

@@ -50,17 +50,18 @@
  *
  *  Calculates the diagonal elements of S as:
  *  \f[
- *  	S_{ii} = factor_ * \sqrt{\frac{4\pi}{a_i}}
+ *  	S_{ii} = factor * \sqrt{\frac{4\pi}{a_i}}
  *  \f]
  *  while the diagonal elements of D are:
  *  \f[
- *  	D_{ii} = -factor_ * \sqrt{\frac{\pi}{a_i}} \frac{1}{R_I}
+ *  	D_{ii} = -factor * \sqrt{\frac{\pi}{a_i}} \frac{1}{R_I}
  *  \f]
  */
 
 struct CollocationIntegrator
 {
-    CollocationIntegrator() : factor_(1.07) {}
+    CollocationIntegrator() : factor(1.07) {}
+    CollocationIntegrator(double f) : factor(f) {}
     ~CollocationIntegrator() {}
 
     /**@{ Single and double layer potentials for a Vacuum Green's function by collocation */
@@ -71,7 +72,7 @@ struct CollocationIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd singleLayer(const Vacuum<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::singleLayer(e,
-                pcm::bind(integrator::SI, this->factor_, 1.0, pcm::_1),
+                pcm::bind(integrator::SI, this->factor, 1.0, pcm::_1),
                 pcm::bind(&Vacuum<DerivativeTraits, CollocationIntegrator>::kernelS, gf, pcm::_1, pcm::_2));
     }
     /*! \tparam DerivativeTraits how the derivatives of the Greens's function are calculated
@@ -81,7 +82,7 @@ struct CollocationIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd doubleLayer(const Vacuum<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::doubleLayer(e,
-                                       pcm::bind(integrator::DI, this->factor_, pcm::_1),
+                                       pcm::bind(integrator::DI, this->factor, pcm::_1),
                                        pcm::bind(&Vacuum<DerivativeTraits, CollocationIntegrator>::kernelD, gf, pcm::_1, pcm::_2, pcm::_3));
     }
     /**@}*/
@@ -94,7 +95,7 @@ struct CollocationIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd singleLayer(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::singleLayer(e,
-                pcm::bind(integrator::SI, this->factor_, gf.epsilon(), pcm::_1),
+                pcm::bind(integrator::SI, this->factor, gf.epsilon(), pcm::_1),
                 pcm::bind(&UniformDielectric<DerivativeTraits, CollocationIntegrator>::kernelS, gf, pcm::_1, pcm::_2));
     }
     /*! \tparam DerivativeTraits how the derivatives of the Greens's function are calculated
@@ -104,7 +105,7 @@ struct CollocationIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd doubleLayer(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::doubleLayer(e,
-                                       pcm::bind(integrator::DI, this->factor_, pcm::_1),
+                                       pcm::bind(integrator::DI, this->factor, pcm::_1),
                                        pcm::bind(&UniformDielectric<DerivativeTraits, CollocationIntegrator>::kernelD, gf, pcm::_1, pcm::_2, pcm::_3));
     }
     /**@}*/
@@ -148,7 +149,7 @@ struct CollocationIntegrator
         for (size_t i = 0; i < mat_size; ++i) {
             // Fill diagonal
             // Diagonal of S inside the cavity
-            double Sii_I = factor_ * std::sqrt(4 * M_PI / e[i].area());
+            double Sii_I = factor * std::sqrt(4 * M_PI / e[i].area());
             // "Diagonal" of Coulomb singularity separation coefficient
             double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
             // "Diagonal" of the image Green's function
@@ -177,9 +178,9 @@ struct CollocationIntegrator
             double area = e[i].area();
             double radius = e[i].sphere().radius;
             // Diagonal of S inside the cavity
-            double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+            double Sii_I = factor * std::sqrt(4 * M_PI / area);
             // Diagonal of D inside the cavity
-            double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
+            double Dii_I = -factor * std::sqrt(M_PI/ area) * (1.0 / radius);
             // "Diagonal" of Coulomb singularity separation coefficient
             double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
             // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
@@ -205,7 +206,7 @@ struct CollocationIntegrator
     /**@}*/
 
     /// Scaling factor for the collocation formulas
-    double factor_;
+    double factor;
 };
 
 #endif // COLLOCATIONINTEGRATOR_HPP

src/bi_operators/NumericalIntegrator.hpp

@@ -53,6 +53,9 @@
 
 struct NumericalIntegrator
 {
+    NumericalIntegrator() : factor(1.07) {}
+    NumericalIntegrator(double f) : factor(f) {}
+    ~NumericalIntegrator() {}
     /**@{ Single and double layer potentials for a Vacuum Green's function by collocation: numerical integration of diagonal */
     /*! \tparam DerivativeTraits how the derivatives of the Greens's function are calculated
      *  \param[in] gf Green's function
@@ -183,6 +186,8 @@ struct NumericalIntegrator
         return Eigen::MatrixXd::Zero(e.size(), e.size());
     }
     /**@}*/
+    /// Scaling factor for the collocation formulas (unused here)
+    double factor;
 };
 
 #endif // NUMERICALINTEGRATOR_HPP

src/bi_operators/PurisimaIntegrator.hpp

@@ -51,7 +51,7 @@
  *
  *  Calculates the diagonal elements of S as:
  *  \f[
- *  	S_{ii} = factor_ * \sqrt{\frac{4\pi}{a_i}}
+ *  	S_{ii} = factor * \sqrt{\frac{4\pi}{a_i}}
  *  \f]
  *  while the diagonal elements of D are:
  *  \f[
@@ -64,7 +64,8 @@
 
 struct PurisimaIntegrator
 {
-    PurisimaIntegrator() : factor_(1.07) {}
+    PurisimaIntegrator() : factor(1.07) {}
+    PurisimaIntegrator(double f) : factor(f) {}
     ~PurisimaIntegrator() {}
 
     /**@{ Single and double layer potentials for a Vacuum Green's function by collocation */
@@ -75,7 +76,7 @@ struct PurisimaIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd singleLayer(const Vacuum<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::singleLayer(e,
-                pcm::bind(integrator::SI, this->factor_, 1.0, pcm::_1),
+                pcm::bind(integrator::SI, this->factor, 1.0, pcm::_1),
                 pcm::bind(&Vacuum<DerivativeTraits, PurisimaIntegrator>::kernelS, gf, pcm::_1, pcm::_2));
     }
     /*! \tparam DerivativeTraits how the derivatives of the Greens's function are calculated
@@ -101,7 +102,7 @@ struct PurisimaIntegrator
     template <typename DerivativeTraits>
     Eigen::MatrixXd singleLayer(const UniformDielectric<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
         return integrator::singleLayer(e,
-                pcm::bind(integrator::SI, this->factor_, gf.epsilon(), pcm::_1),
+                pcm::bind(integrator::SI, this->factor, gf.epsilon(), pcm::_1),
                 pcm::bind(&UniformDielectric<DerivativeTraits, PurisimaIntegrator>::kernelS, gf, pcm::_1, pcm::_2));
     }
     /*! \tparam DerivativeTraits how the derivatives of the Greens's function are calculated
@@ -159,7 +160,7 @@ struct PurisimaIntegrator
         for (size_t i = 0; i < mat_size; ++i) {
             // Fill diagonal
             // Diagonal of S inside the cavity
-            double Sii_I = factor_ * std::sqrt(4 * M_PI / e[i].area());
+            double Sii_I = factor * std::sqrt(4 * M_PI / e[i].area());
             // "Diagonal" of Coulomb singularity separation coefficient
             double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
             // "Diagonal" of the image Green's function
@@ -188,9 +189,9 @@ struct PurisimaIntegrator
             double area = e[i].area();
             double radius = e[i].sphere().radius;
             // Diagonal of S inside the cavity
-            double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+            double Sii_I = factor * std::sqrt(4 * M_PI / area);
             // Diagonal of D inside the cavity
-            double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
+            double Dii_I = -factor * std::sqrt(M_PI/ area) * (1.0 / radius);
             // "Diagonal" of Coulomb singularity separation coefficient
             double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
             // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
@@ -216,7 +217,7 @@ struct PurisimaIntegrator
     /**@}*/
 
     /// Scaling factor for the collocation formulas
-    double factor_;
+    double factor;
     /*! Returns off-diagonal elements of the matrix representation of the double layer operator by collocation
      *  \param[in] elements list of finite elements
      *  \param[in] kernD    function for the evaluation of the off-diagonal of D

src/green/AnisotropicLiquid.hpp

@@ -60,6 +60,12 @@ class AnisotropicLiquid __final : public GreensFunction<DerivativeTraits, Integr
     AnisotropicLiquid(const Eigen::Vector3d & eigen_eps, const Eigen::Vector3d & euler_ang) :
         GreensFunction<DerivativeTraits, IntegratorPolicy, Anisotropic,
                   AnisotropicLiquid<DerivativeTraits, IntegratorPolicy> >() { this->profile_ = Anisotropic(eigen_eps, euler_ang); }
+    /*! \param[in] eigen_eps eigenvalues of the permittivity tensors
+     *  \param[in] euler_ang Euler angles in degrees
+     */
+    AnisotropicLiquid(const Eigen::Vector3d & eigen_eps, const Eigen::Vector3d & euler_ang, double f) :
+        GreensFunction<DerivativeTraits, IntegratorPolicy, Anisotropic,
+                  AnisotropicLiquid<DerivativeTraits, IntegratorPolicy> >(f) { this->profile_ = Anisotropic(eigen_eps, euler_ang); }
     virtual ~AnisotropicLiquid() {}
 
     /*! Calculates the matrix representation of the S operator

src/green/GreenData.hpp

@@ -2,22 +2,22 @@
 /*
  *     PCMSolver, an API for the Polarizable Continuum Model
  *     Copyright (C) 2013-2015 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.readthedocs.org/>
  */
@@ -45,6 +45,8 @@ struct greenData {
     int howProfile;
     /*! The permittivity */
     double epsilon;
+    /*! Scaling for the diagonal of the approximate collocation matrices */
+    double scaling;
     /*! Inverse of the Debye length */
     double kappa;
     /*! Diagonal values of the permittivity tensor with respect to the lab frame */
@@ -77,6 +79,7 @@ struct greenData {
     greenData() { empty = true;}
     greenData(int how_d, int how_i, int how_p,
               double _epsilon = 1.0,
+              double s = 1.07,
               double _kappa = 0.0,
               const Eigen::Vector3d & epstens = Eigen::Vector3d::Zero(),
               const Eigen::Vector3d & euler = Eigen::Vector3d::Zero(),
@@ -91,7 +94,7 @@ struct greenData {
               const Eigen::Vector3d & _o = Eigen::Vector3d::Zero(),
               int l = 30) :
         howDerivative(how_d), howIntegrator(how_i), howProfile(how_p),
-    epsilon(_epsilon), kappa(_kappa), epsilonTensor(epstens), eulerAngles(euler),
+    epsilon(_epsilon), scaling(s), kappa(_kappa), epsilonTensor(epstens), eulerAngles(euler),
 	epsilonReal(_epsReal), epsilonImaginary(_epsImaginary),
     NPspheres(_sphere), NPradii(_sphRadius),
     epsilon1(_e1), epsilon2(_e2), center(_c), width(_w), origin(_o), maxL(l) { empty = false; }

src/green/GreensFunction.hpp

@@ -54,6 +54,7 @@ class GreensFunction: public IGreensFunction
 {
 public:
     GreensFunction() : delta_(1.0e-04), integrator_(IntegratorPolicy()) {}
+    GreensFunction(double f) : delta_(1.0e-04), integrator_(IntegratorPolicy(f)) {}
     virtual ~GreensFunction() {}
     /*! Returns value of the directional derivative of the
      *  Greens's function for the pair of points p1, p2:
@@ -163,6 +164,7 @@ class GreensFunction<Numerical, IntegratorPolicy, ProfilePolicy, Derived>: publi
 {
 public:
     GreensFunction() : delta_(1.0e-04), integrator_(IntegratorPolicy()) {}
+    GreensFunction(double f) : delta_(1.0e-04), integrator_(IntegratorPolicy(f)) {}
     virtual ~GreensFunction() {}
     /*! Returns value of the directional derivative of the
      *  Greens's function for the pair of points p1, p2:

src/green/IonicLiquid.hpp

@@ -57,6 +57,8 @@ class IonicLiquid __final : public GreensFunction<DerivativeTraits, IntegratorPo
 public:
     IonicLiquid(double eps, double k) : GreensFunction<DerivativeTraits, IntegratorPolicy, Yukawa,
                                                   IonicLiquid<DerivativeTraits, IntegratorPolicy> >() { this->profile_ = Yukawa(eps, k); }
+    IonicLiquid(double eps, double k, double f) : GreensFunction<DerivativeTraits, IntegratorPolicy, Yukawa,
+                                                  IonicLiquid<DerivativeTraits, IntegratorPolicy> >(f) { this->profile_ = Yukawa(eps, k); }
     virtual ~IonicLiquid() {}
 
     /*! Calculates the matrix representation of the S operator

src/green/RegisterGreensFunctionToFactory.hpp

@@ -53,8 +53,8 @@ namespace
     struct buildVacuum
     {
         template <typename T, typename U>
-        IGreensFunction * operator()(const greenData & /* data */) {
-            return new Vacuum<T, U>();
+        IGreensFunction * operator()(const greenData & data) {
+            return new Vacuum<T, U>(data.scaling);
         }
     };
 
@@ -74,7 +74,7 @@ namespace
     struct buildUniformDielectric {
         template <typename T, typename U>
         IGreensFunction * operator()(const greenData & data) {
-            return new UniformDielectric<T, U>(data.epsilon);
+            return new UniformDielectric<T, U>(data.epsilon, data.scaling);
         }
     };
 
@@ -94,7 +94,7 @@ namespace
     struct buildIonicLiquid {
         template <typename T, typename U>
         IGreensFunction * operator()(const greenData & data) {
-            return new IonicLiquid<T, U>(data.epsilon, data.kappa);
+            return new IonicLiquid<T, U>(data.epsilon, data.kappa, data.scaling);
         }
     };
 
@@ -114,7 +114,7 @@ namespace
     struct buildAnisotropicLiquid {
         template <typename T, typename U>
         IGreensFunction * operator()(const greenData & data) {
-            return new AnisotropicLiquid<T, U>(data.epsilonTensor, data.eulerAngles);
+            return new AnisotropicLiquid<T, U>(data.epsilonTensor, data.eulerAngles, data.scaling);
         }
     };
 
@@ -134,7 +134,7 @@ namespace
     struct buildSphericalDiffuse {
         template <typename T, typename U>
         IGreensFunction * operator()(const greenData & data) {
-            return new SphericalDiffuse<T, U>(data.epsilon1, data.epsilon2, data.width, data.center, data.origin, data.maxL);
+            return new SphericalDiffuse<T, U>(data.epsilon1, data.epsilon2, data.width, data.center, data.origin, data.maxL, data.scaling);
         }
     };