Solve two linear systems of equations to get ASC when IEFPCM is used

CHANGELOG.md

@@ -6,6 +6,19 @@
 
 - The `CPCMSolver` object now stores the scaled, Hermitian, symmetry-adapted S matrix.
   Polarization weights are then directly computed from the incoming MEP.
+- The `IEFSolver` object now stores the non-Hermitian, symmetry-adapted T and R matrices.
+  The explicit calculation of the inverse of T is thus avoided.
+  However, two square matrices of size equal to the cavity size are stored instead
+  of just one. To obtain the polarization weights _two_ linear systems of equations are solved.
+  A [partially pivoted LU decomposition](http://eigen.tuxfamily.org/dox/classEigen_1_1PartialPivLU.html)
+  is used to solve the linear system(s).
+  The strategy used in v1.1.3 suffered from a reduced numerical accuracy, due to the fact that
+  the polarization weights were not correctly defined.
+
+### Removed
+
+- The `hermitivitize` function will only work correctly on matrices. This
+  reverts modifications in the previous release.
 
 ## [v1.1.3] (2016-07-03)
 

doc/pcmsolver.bib

@@ -299,3 +299,20 @@ @ARTICLE{DiRemigio2016-nn
   doi         = "10.1063/1.4943782",
   addendum        = "Times cited: 0"
 }
+
+@ARTICLE{Barone2004-ae,
+  title   = "{Achieving linear-scaling computational cost for the polarizable
+             continuum model of solvation}",
+  author  = "Barone, Vincenzo and Frisch, Michael J and Scalmani, Giovanni and
+             Kudin, Konstantin N and Scuseria, Gustavo E and Pomelli, Christian
+             S",
+  journal = "Theor. Chem. Acc.",
+  volume  =  111,
+  number  = "2-6",
+  pages   = "90--100",
+  month   =  "1~" # mar,
+  year    =  2004,
+  issn    = "1432-881X, 1432-2234",
+  doi     = "10.1007/s00214-003-0527-2"
+}
+

src/solver/IEFSolver.cpp

@@ -93,8 +93,17 @@ Eigen::VectorXd IEFSolver::computeCharge_impl(const Eigen::VectorXd & potential,
   int irrDim = fullDim/nrBlocks;
   charge.segment(irrep*irrDim, irrDim) =
     - blockTepsilon_[irrep].partialPivLu().solve(blockRinfinity_[irrep] * potential.segment(irrep*irrDim, irrDim));
-  // Symmetrize ASC charge := (charge + charge*)/2
-  if (hermitivitize_) hermitivitize(charge);
+
+  // Obtain polarization weights
+  if (hermitivitize_) {
+    Eigen::VectorXd adj_asc = Eigen::VectorXd::Zero(fullDim);
+    // Form T^\dagger * v = c
+    adj_asc.segment(irrep*irrDim, irrDim) = blockTepsilon_[irrep].adjoint().partialPivLu().solve(potential.segment(irrep*irrDim, irrDim));
+    // Form R^\dagger * c = q^* ("transposed" polarization charges)
+    adj_asc.segment(irrep*irrDim, irrDim) = - blockRinfinity_[irrep].adjoint() * (adj_asc.segment(irrep*irrDim, irrDim).eval());
+    // Get polarization weights
+    charge = 0.5 * (adj_asc + charge.eval());
+  }
 
   return charge;
 }

src/solver/IEFSolver.hpp

@@ -48,7 +48,20 @@ class IGreensFunction;
  *  \note We store the non-Hermitian, symmetry-blocked T(epsilon) and Rinfinity matrices.
  *  The ASC is obtained by multiplying the MEP by Rinfinity and then using a partially
  *  pivoted LU decomposition of T(epsilon) on the resulting vector.
- *  This avoids computing and storing the inverse explicitly.
+ *  In case the polarization weights are requested, we use the approach suggested in
+ *  \cite Barone2004-ae.
+ *  First, the adjoint problem is solved:
+ *  \f[
+ *      \mathbf{T}_\varepsilon^\dagger \tilde{v} = v
+ *  \f]
+ *  Also in this case a partially pivoted LU decomposition is used.
+ *  The "transposed" ASC is obtained by the matrix-vector multiplication:
+ *  \f[
+ *      q^* = \mathbf{R}_\infty^\dagger \tilde{v}
+ *  \f]
+ *  Eventually, the two sets of charges are summed and divided by 2
+ *  This avoids computing and storing the inverse explicitly, at the expense of storing
+ *  both T(epsilon) and Rinfinity.
  */
 
 class IEFSolver : public PCMSolver

src/utils/MathUtils.hpp

@@ -211,17 +211,13 @@ inline void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
  *  \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>
+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.
-  if ((obj_.rows() != 1) && (obj_.cols() != 1)) {
-    obj_ = 0.5 * (obj_ + obj_.adjoint().eval());
-  } else {
-    obj_ = 0.5 * (obj_ + obj_.conjugate().eval());
-  }
+  obj_ = 0.5 * (obj_ + obj_.adjoint().eval());
 }
 
 /*! \fn inline void eulerRotation(Eigen::Matrix3d & R_, const Eigen::Vector3d & eulerAngles_)