Switch to robust LDLt for CPCM.

This avoids computing and storing the explicit inverse of SI.
Symmetrization to polarization weights is done on the computed charges.

Author: Roberto Di Remigio

CHANGELOG.md

@@ -2,6 +2,30 @@
 
 ## [Unreleased]
 
+### Changed
+
+- The `PEDRA.OUT` cavity generator log now reports the initial _and_ final
+  lists of spheres. The final list only contains those spheres that were
+  actually tesselated and, possibly, the added spheres.
+- The `CPCMSolver` object now stores the scaled, unsymmetrized S matrix. The
+  explicit calculation and storage of its inverse is thus avoided. A robust
+  Cholesky decomposition is used to solve the linear equation system. The
+  symmetrization needed to obtain the polarization weights happens directly on
+  the computed charges.
+- The `hermitivitize` function can now correctly symmetrize vectors.
+
+### Fixed
+
+- A fix for the initialization of the explicit list of spheres when running the
+  standalone executable. The bug prevented the generation of the correct
+  `Molecule` object, with subsequent failure in the cavity generator.
+- A memory leak occuring in the cavity generator PEDRA was fixed. This was uncovered by @shoefener
+  and manifested only with considerably large cavities (> 200 input spheres)
+
+### Removed
+
+- The function `CPCMMatrix` in the `SolverImpl.hpp` header file is no longer available.
+
 ## [v1.1.2] (2016-05-31)
 
 ### Fixed
@@ -26,7 +50,7 @@
   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 a feature request
-  from @shofener (Issue #13)
+  from @shoefener (Issue #13)
 - The standalone executable `run_pcm` is now tested within the unit tests
   suite. The tests cover the cases where the cavity is given implicitly,
   explicitly or by substitution of radii on chosen atoms.
@@ -36,12 +60,12 @@
 - 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
+  a feature request from @shoefener (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.
 - CMake minimum required version is 2.8.10
-- Atom, Solvent and Sphere are now PODs. The radii and solvent lists are free
+- `Atom`, `Solvent` and `Sphere` are now PODs. The radii and solvent lists are free
   functions.
 - `PCMSOLVER_ERROR` kills program execution when an error arises but does not
   use C++ exceptions.
@@ -53,15 +77,15 @@
 
 ### Known Issues
 
-- The new printer in GePolCavity might not work properly when an explicit list
+- The new printer in `GePolCavity` might not work properly when an explicit list
   of spheres is provided in the input.
 - On Ubuntu 12.10, 32 bit the Intel compiler version 2013.1 produces a faulty
   library. It is possibly a bug in the implementation of `iso_c_binding`, see
   Issue #25
 
 ### Removed
 
-- SurfaceFunction as a class is no longer available. We keep track of surface
+- `SurfaceFunction` as a class is no longer available. We keep track of surface
   functions at the interface level _via_ a label-vector map.
 
 ## [v1.1.0] (2016-02-07)

doc/users/input.rst

@@ -174,7 +174,7 @@ Medium section keywords
 
 .. glossary::
 
-   Solver
+   SolverType
      Type of solver to be used. All solvers are based on the Integral Equation Formulation of
      the Polarizable Continuum Model :cite:`Cances1998`
 

src/pedra/pedra_cavity.F90

@@ -478,7 +478,7 @@ subroutine polyhedra(intsph,vert,centr,newsph,icav1,icav2,xval,yval, &
                 call output(pts, 1_regint_k, 3_regint_k, 1_regint_k, nv, 3_regint_k, nv, 1_regint_k, lvpri)
             end if
             if(abs(area) <= 1.0d-14) then
-                write(lvpri, '(a, i4)') "Zero area in tessera", nn + 1
+                if (iprpcm >= 10) write(lvpri, '(a, i8)') "Zero area in tessera ", nn + 1
                 go to 310
             end if
             if (nn >= mxts) then
@@ -1313,7 +1313,7 @@ subroutine tessera(ns, nv, pts, ccc, pp, pp1, area, intsph, numts)
         end do
         ICUT = ICUT / 2
         IF(ICUT > 1) THEN
-            WRITE(LVPRI,*) 'Tessera cut in pieces and removed.'
+            IF(IPRPCM >= 10) WRITE(LVPRI,*) 'Tessera cut in pieces and removed.'
             RETURN
         end if
 

src/pedra/pedra_cavity_interface.F90

@@ -166,6 +166,7 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_,         &
 deallocate(centr)
 
 write(pedra_unit, *) "Error code is ", error_code
+write(pedra_unit, *) '<<< Done with PEDRA Fortran code >>>'
 
 close(pedra_unit)
 

src/solver/CPCMSolver.cpp

@@ -32,7 +32,7 @@
 #include "Config.hpp"
 
 #include <Eigen/Core>
-#include <Eigen/LU>
+#include <Eigen/Cholesky>
 
 #include "cavity/Cavity.hpp"
 #include "cavity/Element.hpp"
@@ -43,17 +43,26 @@
 void CPCMSolver::buildSystemMatrix_impl(const Cavity & cavity, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
   if (!isotropic_) PCMSOLVER_ERROR("C-PCM is defined only for isotropic environments!", BOOST_CURRENT_FUNCTION);
-  fullPCMMatrix_ = CPCMMatrix(cavity, gf_i, profiles::epsilon(gf_o.permittivity()), correction_);
-  // Symmetrize K := (K + K+)/2
-  if (hermitivitize_) {
-    hermitivitize(fullPCMMatrix_);
-  }
+  TIMER_ON("Computing SI");
+  double epsilon = profiles::epsilon(gf_o.permittivity());
+  SI_ = gf_i.singleLayer(cavity.elements());
+  SI_ /= (epsilon - 1.0)/(epsilon + correction_);
+  TIMER_OFF("Computing SI");
+
   // Symmetry-pack
   // The number of irreps in the group
   int nrBlocks = cavity.pointGroup().nrIrrep();
   // The size of the irreducible portion of the cavity
   int dimBlock = cavity.irreducible_size();
-  symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
+  // Perform symmetry blocking
+  // If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
+  // into "block diagonal" when all other manipulations are done.
+  if (cavity.pointGroup().nrGenerators() != 0) {
+    TIMER_ON("Symmetry blocking");
+    symmetryBlocking(SI_, cavity.size(), dimBlock, nrBlocks);
+    TIMER_OFF("Symmetry blocking");
+  }
+  symmetryPacking(blockSI_, SI_, dimBlock, nrBlocks);
 
   built_ = true;
 }
@@ -63,12 +72,15 @@ Eigen::VectorXd CPCMSolver::computeCharge_impl(const Eigen::VectorXd & potential
   // The potential and charge vector are of dimension equal to the
   // full dimension of the cavity. We have to select just the part
   // relative to the irrep needed.
-  int fullDim = fullPCMMatrix_.rows();
+  int fullDim = SI_.rows();
   Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
-  int nrBlocks = blockPCMMatrix_.size();
+  int nrBlocks = blockSI_.size();
   int irrDim = fullDim/nrBlocks;
   charge.segment(irrep*irrDim, irrDim) =
-    - blockPCMMatrix_[irrep] * potential.segment(irrep*irrDim, irrDim);
+    - blockSI_[irrep].ldlt().solve(potential.segment(irrep*irrDim, irrDim));
+  // Symmetrize ASC charge := (charge + charge*)/2
+  if (hermitivitize_) hermitivitize(charge);
+
   return charge;
 }
 

src/solver/CPCMSolver.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/>
  */
@@ -43,7 +43,13 @@ class IGreensFunction;
  *  \class CPCMSolver
  *  \brief Solver for conductor-like approximation: C-PCM (COSMO)
  *  \author Roberto Di Remigio
- *  \date 2013
+ *  \date 2013, 2016
+ *
+ *  \note We store the unsymmetrized SI matrix and use a robust
+ *  Cholesky decomposition to solve for the ASC. The ASC is then
+ *  symmetrized. This avoids computing and storing the inverse explicitly.
+ *  The SI matrix is already scaled by the dielectric factor entering the
+ *  definition of the conductor model!
  */
 
 class CPCMSolver : public PCMSolver
@@ -65,10 +71,10 @@ class CPCMSolver : public PCMSolver
     bool hermitivitize_;
     /*! Correction for the conductor results */
     double correction_;
-    /*! PCM matrix, not symmetry blocked */
-    Eigen::MatrixXd fullPCMMatrix_;
-    /*! PCM matrix, symmetry blocked form */
-    std::vector<Eigen::MatrixXd> blockPCMMatrix_;
+    /*! SI matrix, not symmetry blocked */
+    Eigen::MatrixXd SI_;
+    /*! SI matrix, symmetry blocked form */
+    std::vector<Eigen::MatrixXd> blockSI_;
 
     /*! \brief Calculation of the PCM matrix
      *  \param[in] cavity the cavity to be used

src/solver/SolverImpl.hpp

@@ -192,49 +192,6 @@ inline Eigen::MatrixXd isotropicIEFMatrix(const Cavity & cav, const IGreensFunct
   return fullPCMMatrix;
 }
 
-/*! \brief Builds the CPCM matrix
- *  \param[in] cav the discretized cavity
- *  \param[in] gf_i Green's function inside the cavity
- *  \param[in] epsilon permittivity outside the cavity
- *  \param[in] correction CPCM correction factor
- *  \return the \f$ \mathbf{K} = f(\varepsilon)\mathbf{S}^{-1} \f$ matrix
- *
- *  This function calculates the PCM matrix.
- *  The matrix is not symmetrized and is not symmetry packed.
- */
-inline Eigen::MatrixXd CPCMMatrix(const Cavity & cav, const IGreensFunction & gf_i, double epsilon, double correction)
-{
-  // The total size of the cavity
-  PCMSolverIndex cavitySize = cav.size();
-  // The number of irreps in the group
-  int nrBlocks = cav.pointGroup().nrIrrep();
-  // The size of the irreducible portion of the cavity
-  int dimBlock = cav.irreducible_size();
-
-  // Compute SI and DI on the whole cavity, regardless of symmetry
-  TIMER_ON("Computing SI");
-  Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
-  TIMER_OFF("Computing SI");
-
-  // Perform symmetry blocking
-  // If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
-  // into "block diagonal" when all other manipulations are done.
-  if (cav.pointGroup().nrGenerators() != 0) {
-    TIMER_ON("Symmetry blocking");
-    symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-    TIMER_OFF("Symmetry blocking");
-  }
-
-  double fact = (epsilon - 1.0)/(epsilon + correction);
-  // Invert SI  using LU decomposition with full pivoting
-  // This is a rank-revealing LU decomposition, this allows us
-  // to test if SI is invertible before attempting to invert it.
-  Eigen::FullPivLU<Eigen::MatrixXd> SI_LU(SI);
-  if (!(SI_LU.isInvertible()))
-    PCMSOLVER_ERROR("SI matrix is not invertible!", BOOST_CURRENT_FUNCTION);
-  return fact * SI_LU.inverse();
-}
-
 /*! \brief Builds the **anisotropic** \f$ \mathbf{T}_\varepsilon \f$ matrix
  *  \param[in] cav the discretized cavity
  *  \param[in] gf_i Green's function inside the cavity

src/utils/MathUtils.hpp

@@ -203,19 +203,25 @@ inline void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
     }
 }
 
-/*! \fn inline void hermitivitize(Eigen::MatrixBase<Derived> & matrix_)
- *  \param[out] matrix_ the matrix to be hermitivitized
- *  \tparam     Derived the numeric type of matrix_ elements
+/*! \fn inline void hermitivitize(Eigen::MatrixBase<Derived> & obj_)
+ *  \param[out] obj_ the Eigen object to be hermitivitized
+ *  \tparam     Derived the numeric type of obj_ elements
  *
- *  Given matrix_ returns 0.5 * (matrix_ + matrix_^dagger)
+ *  Given obj_ returns 0.5 * (obj_ + obj_^dagger)
+ *  \note We check if a matrix or vector was given, since in the latter
+ *  case we only want the complex conjugation operation to happen.
  */
-template <typename Derived>
-inline void hermitivitize(Eigen::MatrixBase<Derived> & matrix_)
+  template <typename Derived>
+inline void hermitivitize(Eigen::MatrixBase<Derived> & obj_)
 {
-    // We need to use adjoint().eval() to avoid aliasing issues, see:
-    // http://eigen.tuxfamily.org/dox/group__TopicAliasing.html
-    // The adjoint is evaluated explicitly into an intermediate.
-    matrix_ = 0.5 * (matrix_ + matrix_.adjoint().eval());
+  // We need to use adjoint().eval() to avoid aliasing issues, see:
+  // http://eigen.tuxfamily.org/dox/group__TopicAliasing.html
+  // The adjoint is evaluated explicitly into an intermediate.
+  if ((obj_.rows() != 1) && (obj_.cols() != 1)) {
+    obj_ = 0.5 * (obj_ + obj_.adjoint().eval());
+  } else {
+    obj_ = 0.5 * (obj_ + obj_.conjugate().eval());
+  }
 }
 
 /*! \fn inline void eulerRotation(Eigen::Matrix3d & R_, const Eigen::Vector3d & eulerAngles_)