Refactor creation of PCM matrices in the various solvers

Roberto Di Remigio · Roberto Di Remigio · commit c14a8d453bd5 · 2015-11-06T12:57:23.000-05:00
diff --git a/src/solver/CPCMSolver.cpp b/src/solver/CPCMSolver.cpp
@@ -38,69 +38,48 @@
 #include "Element.hpp"
 #include "IGreensFunction.hpp"
 #include "MathUtils.hpp"
+#include "SolverImpl.hpp"
 
 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!");
+  if (!isotropic_) PCMSOLVER_ERROR("C-PCM is defined only for isotropic environments!");
+  fullPCMMatrix_ = CPCMMatrix(cavity, gf_i, profiles::epsilon(gf_o.permittivity()), correction_);
+  // Symmetrize K := (K + K+)/2
+  if (hermitivitize_) {
+    hermitivitize(fullPCMMatrix_);
+  }
+  // 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);
 
-    // The total size of the cavity
-    size_t cavitySize = cavity.size();
-    // 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();
-
-    // Compute SI on the whole cavity, regardless of symmetry
-    Eigen::MatrixXd SI = gf_i.singleLayer(cavity.elements());
-
-    // Perform symmetry blocking only for the SI matrix as the DI matrix is not used.
-    // 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) {
-        symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-    }
-
-    double epsilon = profiles::epsilon(gf_o.permittivity());
-    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!");
-    fullPCMMatrix_ = fact * SI_LU.inverse();
-    // 5. Symmetrize K := (K + K+)/2
-    if (hermitivitize_) {
-        hermitivitize(fullPCMMatrix_);
-    }
-    symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
-
-    built_ = true;
+  built_ = true;
 }
 
 Eigen::VectorXd CPCMSolver::computeCharge_impl(const Eigen::VectorXd & potential, int irrep) const
 {
-    // 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();
-    Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
-    int nrBlocks = blockPCMMatrix_.size();
-    int irrDim = fullDim/nrBlocks;
-    charge.segment(irrep*irrDim, irrDim) =
-        - blockPCMMatrix_[irrep] * potential.segment(irrep*irrDim, irrDim);
-    return charge;
+  // 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();
+  Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
+  int nrBlocks = blockPCMMatrix_.size();
+  int irrDim = fullDim/nrBlocks;
+  charge.segment(irrep*irrDim, irrDim) =
+    - blockPCMMatrix_[irrep] * potential.segment(irrep*irrDim, irrDim);
+  return charge;
 }
 
 std::ostream & CPCMSolver::printSolver(std::ostream & os)
 {
-    os << "Solver Type: C-PCM" << std::endl;
-    if (hermitivitize_) {
-        os << "PCM matrix hermitivitized";
-    } else {
-        os << "PCM matrix NOT hermitivitized (matches old DALTON)";
-    }
+  os << "Solver Type: C-PCM" << std::endl;
+  if (hermitivitize_) {
+    os << "PCM matrix hermitivitized";
+  } else {
+    os << "PCM matrix NOT hermitivitized (matches old DALTON)";
+  }
 
-    return os;
+  return os;
 }
-
diff --git a/src/solver/IEFSolver.cpp b/src/solver/IEFSolver.cpp
@@ -40,6 +40,7 @@
 #include "Element.hpp"
 #include "IGreensFunction.hpp"
 #include "MathUtils.hpp"
+#include "SolverImpl.hpp"
 
 void IEFSolver::buildSystemMatrix_impl(const Cavity & cavity, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
@@ -49,119 +50,36 @@ void IEFSolver::buildSystemMatrix_impl(const Cavity & cavity, const IGreensFunct
 
 void IEFSolver::buildAnisotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
-    // The total size of the cavity
-    size_t 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, DI and SE, DE on the whole cavity, regardless of symmetry
-    Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
-    Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
-    Eigen::MatrixXd SE = gf_o.singleLayer(cav.elements());
-    Eigen::MatrixXd DE = gf_o.doubleLayer(cav.elements());
-
-    // 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) {
-        symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
-        symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-        symmetryBlocking(DE, cavitySize, dimBlock, nrBlocks);
-        symmetryBlocking(SE, cavitySize, dimBlock, nrBlocks);
-    }
-
-    Eigen::MatrixXd a = cav.elementArea().asDiagonal();
-    Eigen::MatrixXd aInv = a.inverse();
-
-    // 1. Form T
-    fullPCMMatrix_ = ((2 * M_PI * aInv - DE) * a * SI + SE * a * (2 * M_PI * aInv + DI.adjoint().eval()));
-    // 2. Invert T using LU decomposition with full pivoting
-    //    This is a rank-revealing LU decomposition, this allows us
-    //    to test if T is invertible before attempting to invert it.
-    Eigen::FullPivLU<Eigen::MatrixXd> T_LU(fullPCMMatrix_);
-    if (!(T_LU.isInvertible())) PCMSOLVER_ERROR("T matrix is not invertible!");
-    fullPCMMatrix_ = T_LU.inverse();
-    Eigen::FullPivLU<Eigen::MatrixXd> SI_LU(SI);
-    if (!(SI_LU.isInvertible())) PCMSOLVER_ERROR("SI matrix is not invertible!");
-    fullPCMMatrix_ *= ((2 * M_PI * aInv - DE) - SE * SI_LU.inverse() * (2 * M_PI * aInv - DI));
-    fullPCMMatrix_ *= a;
-    // 5. Symmetrize K := (K + K+)/2
+    fullPCMMatrix_ = anisotropicIEFMatrix(cav, gf_i, gf_o);
+    // Symmetrize K := (K + K+)/2
     if (hermitivitize_) {
         hermitivitize(fullPCMMatrix_);
     }
     // Pack into a block diagonal matrix
+    // 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();
     // For the moment just packs into a std::vector<Eigen::MatrixXd>
     symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
-    std::ofstream matrixOut("PCM_matrix");
-    matrixOut << "fullPCMMatrix" << std::endl;
-    matrixOut << fullPCMMatrix_ << std::endl;
-    for (int i = 0; i < nrBlocks; ++i) {
-        matrixOut << "Block number " << i << std::endl;
-        matrixOut << blockPCMMatrix_[i] << std::endl;
-    }
 
     built_ = true;
 }
 
 void IEFSolver::buildIsotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
-    // The total size of the cavity
-    size_t 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
-    Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
-    Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
-
-    // 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) {
-        symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
-        symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-    }
-
-    Eigen::MatrixXd a = cav.elementArea().asDiagonal();
-    Eigen::MatrixXd aInv = Eigen::MatrixXd::Zero(cavitySize, cavitySize);
-    aInv = a.inverse();
-
-    // Tq = -Rv -> q = -(T^-1 * R)v = -Kv
-    // T = (2 * M_PI * fact * aInv - DI) * a * SI; R = (2 * M_PI * aInv - DI)
-    // fullPCMMatrix_ = K = T^-1 * R * a
-    // 1. Form T
-    double epsilon = profiles::epsilon(gf_o.permittivity());
-    double fact = (epsilon + 1.0)/(epsilon - 1.0);
-    fullPCMMatrix_ = (2 * M_PI * fact * aInv - DI) * a * SI;
-    // 2. Invert T using LU decomposition with full pivoting
-    //    This is a rank-revealing LU decomposition, this allows us
-    //    to test if T is invertible before attempting to invert it.
-    Eigen::FullPivLU<Eigen::MatrixXd> T_LU(fullPCMMatrix_);
-    if (!(T_LU.isInvertible()))
-        PCMSOLVER_ERROR("T matrix is not invertible!");
-    fullPCMMatrix_ = T_LU.inverse();
-    // 3. Multiply T^-1 and R
-    fullPCMMatrix_ *= (2 * M_PI * aInv - DI);
-    // 4. Multiply by a
-    fullPCMMatrix_ *= a;
-    // 5. Symmetrize K := (K + K+)/2
+    fullPCMMatrix_ = isotropicIEFMatrix(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
+    // Symmetrize K := (K + K+)/2
     if (hermitivitize_) {
         hermitivitize(fullPCMMatrix_);
     }
     // Pack into a block diagonal matrix
+    // 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();
     // For the moment just packs into a std::vector<Eigen::MatrixXd>
     symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
-    std::ofstream matrixOut("PCM_matrix");
-    matrixOut << "fullPCMMatrix" << std::endl;
-    matrixOut << fullPCMMatrix_ << std::endl;
-    for (int i = 0; i < nrBlocks; ++i) {
-        matrixOut << "Block number " << i << std::endl;
-        matrixOut << blockPCMMatrix_[i] << std::endl;
-    }
 
     built_ = true;
 }
diff --git a/src/solver/PCMSolver.hpp b/src/solver/PCMSolver.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/>
  */
diff --git a/src/solver/SolverImpl.hpp b/src/solver/SolverImpl.hpp
