Fix bug and improve ODE integration of the spherical diffuse Green's …

…function (#109)

The integration of the differential equation for the angular momentum component is now performed after a second change of variables. The asymptotic solution is no longer logarithmic but linear.
The component at L=0 has now been reintroduced. Without it, we don't get the correct solvation energy.

.gitignore

@@ -37,3 +37,7 @@
 
 # Autocmake stuff
 cmake/lib/*.pyc
+
+# Emacs backup files
+*~
+

CHANGELOG.md

@@ -4,6 +4,19 @@
 
 ### Added
 
+- Double logarithmic scale for the integration of spherical diffuse
+  interfaces: much more stable than the previous version, allowing for
+  Runge-Kutta 4 integrator.
+
+### Fixed
+
+- Bug in the diffuse interface Green's function. Contrary to the sharp
+  interface case, it is wrong to remove the monopole, which becomes
+  identically zero when the corresponding differential equation is
+  solved in extreme cases (e.g. charge far away from the sphere).
+
+### Added
+
 - A new CMake module `options_wrappers.cmake` that adds new wrapper macros for
   the CMake `option` command.
 
@@ -125,6 +138,8 @@
 - The uppercased contents of the `.pcm` input file are written to a temporary
   file, instead of overwriting the user provided file. The temporary file is
   removed after it has been parsed. Fixes #91 as noted by @ilfreddy.
+- Use Runge-Kutta-Fehlberg 7(8) ODE solver to integrate the radial equation
+  in the spherical diffuse Green's function class.
 
 ### Fixed
 

doc/code-reference/greens-functions.rst

@@ -7,19 +7,12 @@ Green's functions objects relies on the Factory Method pattern
 :cite:`Gamma1994,Alexandrescu2001`, implemented through the
 generic Factory class.
 
-The top-level header, _i.e._ to be included in client code, is `Green.hpp`.
-The common interface to all Green's function classes is specified by the `IGreensFunction` class,
+The top-level header, _i.e._ to be included in client code, is ``Green.hpp``.
+The common interface to all Green's function classes is specified by the ``IGreensFunction`` class,
 this is non-templated.
-All other classes are templated. TODO: explain template parameters.
+All other classes are templated.
 The Green's functions are registered to the factory based on a label encoding: type, derivative, and dielectric profile.
-The only allowed labels are:
-- `VACUUM_NUMERICAL`
-- `VACUUM_DERIVATIVE`
-- `VACUMM_GRADIENT`
-- `VACUUM_HESSIAN`
-
-- `SPHERICALDIFFUSE_TANH`
-- `SPHERICALDIFFUSE_ERF`
+The only allowed labels must be listed in ``src/green/Green.hpp``. If they are not, they can not be selected at run time.
 
 .. image:: ../gfx/green.png
    :scale: 70 %

doc/pcmsolver.bib

@@ -329,3 +329,13 @@ @ARTICLE{Allinger1994-tx
   doi     = "10.1016/S0166-1280(09)80008-0"
 }
 
+@article{Fosso-Tande2013,
+title = "Implicit solvation models in a multiresolution multiwavelet basis",
+journal = "Chemical Physics Letters",
+volume = "561-562",
+pages = "179--184",
+year = "2013",
+issn = "0009-2614",
+doi = "https://doi.org/10.1016/j.cplett.2013.01.065",
+author = "Jacob Fosso-Tande and Robert J. Harrison"
+}

doc/users/input.rst

@@ -363,13 +363,13 @@ while the Green's function outside might vary.
      * **Valid values**: :math:`\varepsilon \geq 1.0`
      * **Default**: 1.0
 
-    Profile
+   Profile
       Functional form of the dielectric profile
 
       * **Type**: string
-      * **Valid values**: Tanh | Erf
+      * **Valid values**: Tanh | Erf | Log
       * **Valid for**: SphericalDiffuse
-      * **Default**: Tanh
+      * **Default**: Log
 
    Eps1
      Static dielectric permittivity inside the interface

src/green/Green.hpp

@@ -33,7 +33,6 @@
 #include "SphericalDiffuse.hpp"
 #include "UniformDielectric.hpp"
 #include "Vacuum.hpp"
-#include "dielectric_profile/MembraneTanh.hpp"
 #include "dielectric_profile/OneLayerErf.hpp"
 #include "dielectric_profile/OneLayerTanh.hpp"
 #include "utils/Factory.hpp"
@@ -89,6 +88,8 @@ inline Factory<detail::CreateGreensFunction> bootstrapFactory() {
                      createSphericalDiffuse<dielectric_profile::OneLayerTanh>);
   factory_.subscribe("SPHERICALDIFFUSE_NUMERICAL_ERF",
                      createSphericalDiffuse<dielectric_profile::OneLayerErf>);
+  factory_.subscribe("SPHERICALDIFFUSE_NUMERICAL_LOG",
+                     createSphericalDiffuse<dielectric_profile::OneLayerLog>);
 
   return factory_;
 }

src/green/GreensFunction.hpp

@@ -149,6 +149,7 @@ class GreensFunction : public IGreensFunction {
    */
   virtual DerivativeTraits operator()(DerivativeTraits * source,
                                       DerivativeTraits * probe) const = 0;
+
   /*! Returns value of the kernel of the \f$\mathcal{S}\f$ integral operator, i.e.
    *  the value of the Greens's function for the pair of points p1, p2:
    *  \f$ G(\mathbf{p}_1, \mathbf{p}_2)\f$
@@ -169,6 +170,7 @@ class GreensFunction : public IGreensFunction {
     pp[2] = p2(2);
     return this->operator()(sp, pp)[0];
   }
+
   virtual std::ostream & printObject(std::ostream & os) __override {
     os << "Green's Function" << std::endl;
     return os;