Merge remote-tracking branch 'upstream/release/1.Y' into release/1.Y

Author: robertodr

.gitattributes

@@ -67,6 +67,7 @@ tests/green/*.cpp                licensefile=.githooks/LICENSE-C++
 tests/iefpcm/*.cpp               licensefile=.githooks/LICENSE-C++
 tests/input/*.cpp                licensefile=.githooks/LICENSE-C++
 tests/numerical_quadrature/*.cpp licensefile=.githooks/LICENSE-C++
+tests/utils/*.cpp                licensefile=.githooks/LICENSE-C++
 
 # tools
 doc/conf.py               licensefile=.githooks/LICENSE-Python

CHANGELOG.md

@@ -7,6 +7,10 @@
 - A Python script, using Matplotlib, to plot the cavity.
   The script can also color-map the finite elements according to the values of
   a surface function.
+- The input learnt to parse the additional `ChargeDistribution` section.
+  It is possible to specify a classical charge distribution of point multipoles.
+  This can be an additional source of electrostatic potential for the calculation
+  of the ASC.
 
 ## [Version 1.1.10] - 2017-03-27
 

api/pcmsolver.f90

@@ -187,7 +187,7 @@ end subroutine pcmsolver_set_surface_function
     subroutine pcmsolver_print_surface_function(context, name) bind(C)
       import
       type(c_ptr), value :: context
-      character(c_char), intent(in) :: name
+      character(c_char), intent(in) :: name(*)
     end subroutine pcmsolver_print_surface_function
 
     subroutine pcmsolver_save_surface_functions(context) bind(C)
@@ -198,13 +198,13 @@ end subroutine pcmsolver_save_surface_functions
     subroutine pcmsolver_save_surface_function(context, name) bind(C)
       import
       type(c_ptr), value :: context
-      character(kind=c_char, len=1), intent(in) :: name
+      character(kind=c_char, len=1), intent(in) :: name(*)
     end subroutine pcmsolver_save_surface_function
 
     subroutine pcmsolver_load_surface_function(context, name) bind(C)
       import
       type(c_ptr), value :: context
-      character(kind=c_char, len=1), intent(in) :: name
+      character(kind=c_char, len=1), intent(in) :: name(*)
     end subroutine pcmsolver_load_surface_function
 
     subroutine pcmsolver_write_timings(context) bind(C)

api/pcmsolver.h

@@ -246,13 +246,19 @@ PCMSolver_API void pcmsolver_save_surface_functions(pcmsolver_context_t * contex
 /*! \brief Dumps a surface function to NumPy array in .npy file
  *  \param[in, out] context the PCM context object
  *  \param[in] name label of the surface function
+ *
+ *  \note The name parameter is the name of the NumPy array file
+ *  **without** .npy extension
  */
 PCMSolver_API void pcmsolver_save_surface_function(pcmsolver_context_t * context,
                                                    const char * name);
 
 /*! \brief Loads a surface function from a .npy file
  *  \param[in, out] context the PCM context object
  *  \param[in] name label of the surface function
+ *
+ *  \note The name parameter is the name of the NumPy array file
+ *  **without** .npy extension
  */
 PCMSolver_API void pcmsolver_load_surface_function(pcmsolver_context_t * context,
                                                    const char * name);

doc/code-reference/helper-classes.rst

@@ -29,6 +29,14 @@ Atom
    :protected-members:
    :private-members:
 
+ChargeDistribution
+------------------
+.. doxygenstruct:: ChargeDistribution
+   :project: PCMSolver
+   :members:
+   :protected-members:
+   :private-members:
+
 Molecule
 --------
 .. doxygenclass:: Molecule

doc/snippets/example_input-all.inp

@@ -13,10 +13,13 @@ Cavity {
        Spheres = [Array of Doubles]
 }
 Medium {
+       Nonequilibrium = [Bool]
        Solvent = [String]
        SolverType = [String]
        MatrixSymm = [Bool]
        Correction = [Double]
+       DiagonalIntegrator = [String]
+       DiagonalScaling = [Double]
        ProbeRadius = [Double]
        Green<GreenTag> {
              Type = [String]
@@ -33,3 +36,11 @@ Medium {
              MaxL = [Integer]
        }
 }
+Molecule {
+   MEP = [Bool]
+   Geometry = [Double]
+}
+ChargeDistribution {
+  Monopoles = [Double]
+  Dipoles = [Double]
+}

doc/users/input.rst

@@ -15,6 +15,21 @@ core parameters.
 In this page, input style and input parameters available in Method 1 will be
 documented.
 
+Note that it is also possible to run the module standalone and use a classical charge
+distribution.
+The classical charge distribution can be specified by giving a molecular geometry
+in the molecule section and an additional point multipoles distribution
+in the charge distribution section.
+The ``run_pcm.x`` executable has to be compiled for a standalone run with:
+
+.. code-block:: bash
+
+   python pcmsolver.py -x molecule.inp
+
+where the ``molecule.inp`` input file looks like:
+
+.. literalinclude:: ../snippets/molecule.inp
+
 Input style
 -----------
 
@@ -46,6 +61,14 @@ Available sections:
   i.e. the Green's functions inside and outside the cavity;
 + Green, subsection of medium. Sets up the Green's function inside and
   outside the cavity.
++ Molecule: molecular geometry to be used in a standalone run.
++ ChargeDistribution: sets up a classical multipolar (currently up to dipoles)
+  charge distribution to use as additional source of electrostatic potential.
+
+.. note::
+
+   The Molecule and ChargeDistribution sections only make sense in a standalone run,
+   i.e. when using the ``run_pcm.x`` executable.
 
 .. warning::
 
@@ -436,6 +459,27 @@ where the ``molecule.inp`` input file looks like:
 
      * **Type**: array of doubles
 
+ChargeDistribution section keywords
+-----------------------------------
+
+Set a classical charge distribution, inside or outside the cavity
+No additional spheres will be generated.
+
+.. glossary::
+
+   Monopoles
+     Array of point charges
+     Format is :math:`[\ldots, x_i, y_i, z_i, Q_i, \ldots]`
+
+     * **Type**: array of doubles
+
+   Dipoles
+     Array of point dipoles.
+     Format is :math:`[\ldots, x_i, y_i, z_i, \mu_{x_i}, \mu_{y_i}, \mu_{z_i} \ldots]`
+     The dipole moment components are always read in atomic units.
+
+     * **Type**: array of doubles
+
 .. _available-radii:
 
 Available radii

src/bin/run_pcm.cpp

@@ -41,6 +41,7 @@
 #include "interface/Input.hpp"
 #include "interface/Meddle.hpp"
 #include "utils/Molecule.hpp"
+#include "utils/ChargeDistribution.hpp"
 
 std::ofstream pcmsolver_out;
 
@@ -66,16 +67,17 @@ int main(int argc, char * argv[]) {
   PCMSolverIndex size = context_.getCavitySize();
 
   // Form vector with electrostatic potential
-  // FIXME
-  // 1. Try to understand why this is needed
-  // 2. Try to re-write this such that the input object is not needed!!!
   // First compute the potential from the classical point multipoles distribution
   // then add the one from the molecule
   TIMER_ON("Computing MEP");
-  // Eigen::VectorXd mep = Eigen::VectorXd::Zero(size);
-  // if (input.MEPfromMolecule())
-  //  mep += computeMEP(context_.molecule(), context_.getCenters());
-  Eigen::VectorXd mep = computeMEP(context_.molecule(), context_.getCenters());
+  // FIXME currently hardcoded to the dipole-dipole interaction potential in vacuum
+  Eigen::VectorXd mep =
+      computeDipolarPotential(context_.getCenters(), input.multipoles());
+  // FIXME
+  // 1. Try to understand why this is needed
+  // 2. Try to re-write this such that the input object is not needed!!!
+  if (input.MEPfromMolecule())
+    mep += computeMEP(context_.molecule(), context_.getCenters());
   TIMER_OFF("Computing MEP");
   context_.setSurfaceFunction(mep.size(), mep.data(), "MEP");
   // Compute apparent surface charge

src/green/AnisotropicLiquid.cpp

@@ -93,6 +93,16 @@ KernelD AnisotropicLiquid<DerivativeTraits>::exportKernelD_impl() const {
                    pcm::_3);
 }
 
+template <typename DerivativeTraits>
+DerivativeProbe AnisotropicLiquid<DerivativeTraits>::exportDerivativeProbe_impl()
+    const {
+  return pcm::bind(&AnisotropicLiquid<DerivativeTraits>::derivativeProbe,
+                   *this,
+                   pcm::_1,
+                   pcm::_2,
+                   pcm::_3);
+}
+
 template <typename DerivativeTraits>
 double AnisotropicLiquid<DerivativeTraits>::singleLayer_impl(
     const Element & /* e */,

src/green/AnisotropicLiquid.hpp

@@ -78,6 +78,7 @@ class AnisotropicLiquid __final
 
   virtual KernelS exportKernelS_impl() const __override;
   virtual KernelD exportKernelD_impl() const __override;
+  virtual DerivativeProbe exportDerivativeProbe_impl() const __override;
 
   __noreturn virtual double singleLayer_impl(const Element & /* e */,
                                              double /* factor */) const __override;