Expose proper interface for MMFQ

Author: robertodr

src/bin/run_pcm.cpp

@@ -70,12 +70,12 @@ int main(int argc, char * argv[]) {
   // First compute the potential from the classical point multipoles distribution
   // then add the one from the molecule
   TIMER_ON("Computing MEP");
-  // 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
+  // 1. Currently hardcoded to the dipole-dipole interaction potential in vacuum
   // 2. Try to re-write this such that the input object is not needed!!!
+  Eigen::VectorXd mep = Eigen::VectorXd::Zero(size);
+  if (input.MEPfromChargeDistribution())
+    mep += computeDipolarPotential(context_.getCenters(), input.multipoles());
   if (input.MEPfromMolecule())
     mep += computeMEP(context_.molecule(), context_.getCenters());
   TIMER_OFF("Computing MEP");
@@ -90,6 +90,7 @@ int main(int argc, char * argv[]) {
   Eigen::VectorXd rsp_asc(size);
   context_.getSurfaceFunction(rsp_asc.size(), rsp_asc.data(), "RspASC");
   TIMER_OFF("Computing ASC");
+  context_.printSurfaceFunction("ASC");
   // Compute energy and print it out
   pcmsolver_out << "Solvation energy = "
                 << std::setprecision(std::numeric_limits<long double>::digits10)

src/interface/Input.cpp

@@ -160,7 +160,9 @@ void Input::reader(const std::string & filename) {
   isDynamic_ = medium.getBool("NONEQUILIBRIUM");
 
   const Section & chgdist = input_.getSect("CHARGEDISTRIBUTION");
+  MEPfromChargeDist_ = false;
   if (chgdist.isDefined()) {
+    MEPfromChargeDist_ = true;
     // Set monopoles
     if (chgdist.getKey<std::vector<double> >("MONOPOLES").isDefined()) {
       std::vector<double> mono = chgdist.getDblVec("MONOPOLES");
@@ -190,26 +192,27 @@ void Input::reader(const std::string & filename) {
         j += 6;
       }
     }
-    // Set fluctuating charges
-    isFQ_ = false;
-    const Section & mmfq = input_.getSect("MMFQ");
-    if (mmfq.isDefined()) {
-      isFQ_ = true;
-      fragments_.nSitesPerFragment =
-          static_cast<PCMSolverIndex>(mmfq.getInt("SITESPERFRAGMENT"));
-      std::vector<double> sites = mmfq.getDblVec("SITES");
-      int j = 0;
-      int n = int(sites.size() / 5);
-      fragments_.sites = Eigen::Matrix3Xd::Zero(3, n);
-      fragments_.chi = Eigen::VectorXd::Zero(n);
-      fragments_.eta = Eigen::VectorXd::Zero(n);
-      for (int i = 0; i < n; ++i) {
-        fragments_.sites.col(i) =
-            (Eigen::Vector3d() << sites[j], sites[j + 1], sites[j + 2]).finished();
-        fragments_.chi(i) = sites[j + 3];
-        fragments_.eta(i) = sites[j + 4];
-        j += 5;
-      }
+  }
+  // Set fluctuating charges
+  isFQ_ = false;
+  const Section & mmfq = input_.getSect("MMFQ");
+  if (mmfq.isDefined()) {
+    isFQ_ = true;
+    fragments_.nSitesPerFragment =
+        static_cast<PCMSolverIndex>(mmfq.getInt("SITESPERFRAGMENT"));
+    std::vector<double> sites = mmfq.getDblVec("SITES");
+    int j = 0;
+    int n = int(sites.size() / 5);
+    fragments_.nFragments = int(n / fragments_.nSitesPerFragment);
+    fragments_.sites = Eigen::Matrix3Xd::Zero(3, n);
+    fragments_.chi = Eigen::VectorXd::Zero(n);
+    fragments_.eta = Eigen::VectorXd::Zero(n);
+    for (int i = 0; i < n; ++i) {
+      fragments_.sites.col(i) =
+          (Eigen::Vector3d() << sites[j], sites[j + 1], sites[j + 2]).finished();
+      fragments_.chi(i) = sites[j + 3];
+      fragments_.eta(i) = sites[j + 4];
+      j += 5;
     }
   }
 

src/interface/Input.hpp

@@ -136,6 +136,7 @@ class PCMSolver_EXPORT Input {
   ChargeDistribution multipoles() const { return multipoles_; }
   MMFQ fragments() const { return fragments_; }
   bool MEPfromMolecule() { return MEPfromMolecule_; }
+  bool MEPfromChargeDistribution() { return MEPfromChargeDist_; }
 
   /// Operators
   /// operator<<
@@ -249,6 +250,8 @@ class PCMSolver_EXPORT Input {
   bool isFQ_;
   /// Whether to calculate the MEP from the molecular geometry
   bool MEPfromMolecule_;
+  /// Whether to calculate the MEP from the charge distribution
+  bool MEPfromChargeDist_;
   /// Classical charge distribution of point multipoles
   ChargeDistribution multipoles_;
   /// Classical fluctuating charges MM force field

src/interface/Meddle.cpp

@@ -145,18 +145,25 @@ void Meddle::CTORBody() {
               << std::endl;
   infoStream_ << "Input parsing done " << input_.providedBy() << std::endl;
 
-  TIMER_ON("Meddle::initCavity");
-  initCavity();
-  TIMER_OFF("Meddle::initCavity");
-
-  TIMER_ON("Meddle::initStaticSolver");
-  initStaticSolver();
-  TIMER_OFF("Meddle::initStaticSolver");
-
-  if (input_.isDynamic()) {
-    TIMER_ON("Meddle::initDynamicSolver");
-    initDynamicSolver();
-    TIMER_OFF("Meddle::initDynamicSolver");
+  if (input_.isFQ()) { /* MMFQ calculation */
+    hasFQ_ = true;
+    TIMER_ON("Meddle::initMMFQ");
+    initMMFQ();
+    TIMER_OFF("Meddle::initMMFQ");
+  } else { /* Pure PCM calculation */
+    TIMER_ON("Meddle::initCavity");
+    initCavity();
+    TIMER_OFF("Meddle::initCavity");
+
+    TIMER_ON("Meddle::initStaticSolver");
+    initStaticSolver();
+    TIMER_OFF("Meddle::initStaticSolver");
+
+    if (input_.isDynamic()) {
+      TIMER_ON("Meddle::initDynamicSolver");
+      initDynamicSolver();
+      TIMER_OFF("Meddle::initDynamicSolver");
+    }
   }
 }
 
@@ -167,7 +174,8 @@ Meddle::Meddle(const Input & input, const HostWriter & write)
       cavity_(nullptr),
       K_0_(nullptr),
       K_d_(nullptr),
-      hasDynamic_(false) {
+      hasDynamic_(false),
+      hasFQ_(false) {
   input_.initMolecule();
   CTORBody();
 }
@@ -179,7 +187,8 @@ Meddle::Meddle(const std::string & inputFileName, const HostWriter & write)
       cavity_(nullptr),
       K_0_(nullptr),
       K_d_(nullptr),
-      hasDynamic_(false) {
+      hasDynamic_(false),
+      hasFQ_(false) {
   input_.initMolecule();
   CTORBody();
 }
@@ -196,7 +205,8 @@ Meddle::Meddle(int nr_nuclei,
       cavity_(nullptr),
       K_0_(nullptr),
       K_d_(nullptr),
-      hasDynamic_(false) {
+      hasDynamic_(false),
+      hasFQ_(false) {
   TIMER_ON("Meddle::initInput");
   initInput(nr_nuclei, charges, coordinates, symmetry_info);
   TIMER_OFF("Meddle::initInput");
@@ -216,7 +226,8 @@ Meddle::Meddle(int nr_nuclei,
       cavity_(nullptr),
       K_0_(nullptr),
       K_d_(nullptr),
-      hasDynamic_(false) {
+      hasDynamic_(false),
+      hasFQ_(false) {
   TIMER_ON("Meddle::initInput");
   initInput(nr_nuclei, charges, coordinates, symmetry_info);
   TIMER_OFF("Meddle::initInput");
@@ -235,7 +246,8 @@ Meddle::Meddle(int nr_nuclei,
       cavity_(nullptr),
       K_0_(nullptr),
       K_d_(nullptr),
-      hasDynamic_(false) {
+      hasDynamic_(false),
+      hasFQ_(false) {
   // This one does a deferred initialization:
   // The CTORBody function is not called at this point, but during refresh
   TIMER_ON("Meddle::initInput");
@@ -250,22 +262,26 @@ void pcmsolver_delete(pcmsolver_context_t * context) {
   delete AS_TYPE(pcm::Meddle, context);
 }
 pcm::Meddle::~Meddle() {
-  delete cavity_;
-  delete K_0_;
-  if (hasDynamic_)
-    delete K_d_;
+  if (hasFQ_) { /* MMFQ calculation */
+    delete FQ_;
+  } else { /* Pure PCM calculation */
+    delete cavity_;
+    delete K_0_;
+    if (hasDynamic_)
+      delete K_d_;
+  }
 }
 
 PCMSolverIndex pcmsolver_get_cavity_size(pcmsolver_context_t * context) {
   return (AS_CTYPE(pcm::Meddle, context)->getCavitySize());
 }
-PCMSolverIndex pcm::Meddle::getCavitySize() const { return cavity_->size(); }
+PCMSolverIndex pcm::Meddle::getCavitySize() const { return std::get<0>(size_); }
 
 PCMSolverIndex pcmsolver_get_irreducible_cavity_size(pcmsolver_context_t * context) {
   return (AS_CTYPE(pcm::Meddle, context)->getIrreducibleCavitySize());
 }
 PCMSolverIndex pcm::Meddle::getIrreducibleCavitySize() const {
-  return cavity_->irreducible_size();
+  return std::get<1>(size_);
 }
 
 void pcmsolver_get_centers(pcmsolver_context_t * context, double centers[]) {
@@ -275,16 +291,26 @@ void pcmsolver_get_centers(pcmsolver_context_t * context, double centers[]) {
 }
 void pcm::Meddle::getCenters(double centers[]) const {
   TIMER_ON("Meddle::getCenters");
-  Eigen::Map<Eigen::Matrix3Xd>(centers, 3, cavity_->size()) =
-      cavity_->elementCenter();
+  if (hasFQ_) { /* MMFQ calculation */
+    PCMSolverIndex size = input_.fragments().sites.cols();
+    Eigen::Map<Eigen::Matrix3Xd>(centers, 3, size) = input_.fragments().sites;
+  } else { /* Pure PCM calculation */
+    Eigen::Map<Eigen::Matrix3Xd>(centers, 3, cavity_->size()) =
+        cavity_->elementCenter();
+  }
   TIMER_OFF("Meddle::getCenters");
 }
 
 void pcmsolver_get_center(pcmsolver_context_t * context, int its, double center[]) {
   AS_CTYPE(pcm::Meddle, context)->getCenter(its, center);
 }
 void pcm::Meddle::getCenter(int its, double center[]) const {
-  Eigen::Map<Eigen::Vector3d>(center, 3, 1) = cavity_->elementCenter(its - 1);
+  if (hasFQ_) { /* MMFQ calculation */
+    Eigen::Map<Eigen::Matrix3Xd>(center, 3, 1) =
+        input_.fragments().sites.col(its - 1);
+  } else { /* Pure PCM calculation */
+    Eigen::Map<Eigen::Vector3d>(center, 3, 1) = cavity_->elementCenter(its - 1);
+  }
 }
 
 void pcmsolver_get_areas(pcmsolver_context_t * context, double areas[]) {
@@ -303,7 +329,16 @@ double pcmsolver_compute_polarization_energy(pcmsolver_context_t * context,
 }
 double pcm::Meddle::computePolarizationEnergy(const std::string & mep_name,
                                               const std::string & asc_name) const {
-  double energy = functions_.at(mep_name).dot(functions_.at(asc_name));
+  double energy = 0.0;
+  if (hasFQ_) { /* MMFQ calculation */
+    // Dot product of MEP + Electronegativities times Fluctuating charges
+    energy = (functions_.at(mep_name) + input_.fragments().chi)
+                 .dot(functions_.at(asc_name));
+  } else { /* Pure PCM calculation */
+    // Dot product of MEP and ASC surface function
+    energy =
+        functions_.at(mep_name).dot(functions_.at(asc_name));
+  }
   return (energy / 2.0);
 }
 
@@ -335,9 +370,14 @@ void pcm::Meddle::computeASC(const std::string & mep_name,
                              int irrep) {
   // Get the proper iterators
   SurfaceFunctionMapConstIter iter_pot = functions_.find(mep_name);
-  Eigen::VectorXd asc = K_0_->computeCharge(iter_pot->second, irrep);
-  // Renormalize
-  asc /= double(cavity_->pointGroup().nrIrrep());
+  Eigen::VectorXd asc = Eigen::VectorXd::Zero(iter_pot->second.size());
+  if (hasFQ_) { /* MMFQ calculation */
+    asc = FQ_->computeCharge(iter_pot->second);
+  } else { /* Pure PCM calculation */
+    asc = K_0_->computeCharge(iter_pot->second, irrep);
+    // Renormalize
+    asc /= double(cavity_->pointGroup().nrIrrep());
+  }
   // Insert it into the map
   if (functions_.count(asc_name) == 1) { // Key in map already
     functions_[asc_name] = asc;
@@ -360,14 +400,18 @@ void pcm::Meddle::computeResponseASC(const std::string & mep_name,
                                      int irrep) {
   // Get the proper iterators
   SurfaceFunctionMapConstIter iter_pot = functions_.find(mep_name);
-  Eigen::VectorXd asc(cavity_->size());
-  if (hasDynamic_) {
-    asc = K_d_->computeCharge(iter_pot->second, irrep);
-  } else {
-    asc = K_0_->computeCharge(iter_pot->second, irrep);
+  Eigen::VectorXd asc = Eigen::VectorXd::Zero(iter_pot->second.size());
+  if (hasFQ_) { /* MMFQ calculation */
+    asc = FQ_->computeCharge(iter_pot->second);
+  } else { /* Pure PCM calculation */
+    if (hasDynamic_) {
+      asc = K_d_->computeCharge(iter_pot->second, irrep);
+    } else {
+      asc = K_0_->computeCharge(iter_pot->second, irrep);
+    }
+    // Renormalize
+    asc /= double(cavity_->pointGroup().nrIrrep());
   }
-  // Renormalize
-  asc /= double(cavity_->pointGroup().nrIrrep());
   if (functions_.count(asc_name) == 1) { // Key in map already
     functions_[asc_name] = asc;
   } else { // Create key-value pair
@@ -387,7 +431,7 @@ void pcmsolver_get_surface_function(pcmsolver_context_t * context,
 void pcm::Meddle::getSurfaceFunction(PCMSolverIndex size,
                                      double values[],
                                      const std::string & name) const {
-  if (cavity_->size() != size)
+  if (std::get<0>(size_) != size)
     PCMSOLVER_ERROR("The " + name + " SurfaceFunction is bigger than the cavity!");
 
   SurfaceFunctionMapConstIter iter = functions_.find(name);
@@ -408,7 +452,7 @@ void pcmsolver_set_surface_function(pcmsolver_context_t * context,
 void pcm::Meddle::setSurfaceFunction(PCMSolverIndex size,
                                      double values[],
                                      const std::string & name) {
-  if (cavity_->size() != size)
+  if (std::get<0>(size_) != size)
     PCMSOLVER_ERROR("The " + name + " SurfaceFunction is bigger than the cavity!");
 
   Eigen::VectorXd func = Eigen::Map<Eigen::VectorXd>(values, size, 1);
@@ -460,7 +504,7 @@ void pcmsolver_load_surface_function(pcmsolver_context_t * context,
 void pcm::Meddle::loadSurfaceFunction(const std::string & name) {
   hostWriter_("\nLoading surface function " + name + " from .npy file");
   Eigen::VectorXd values = cnpy::custom::npy_load<double>(name + ".npy");
-  if (values.size() != cavity_->size())
+  if (values.size() != std::get<0>(size_))
     PCMSOLVER_ERROR("The loaded " + name +
                     " surface function is bigger than the cavity!");
   // Append to global map
@@ -656,6 +700,7 @@ void Meddle::initInput(int nr_nuclei,
 void Meddle::initCavity() {
   cavity_ = cavity::bootstrapFactory().create(input_.cavityParams().cavityType,
                                               input_.cavityParams());
+  size_ = std::make_tuple(cavity_->size(), cavity_->irreducible_size());
   cavity_->saveCavity();
 
   infoStream_ << "========== Cavity " << std::endl;
@@ -679,7 +724,7 @@ void Meddle::initStaticSolver() {
   delete biop;
 
   // Perform Gauss' theorem check for nuclear charges
-  GaussCheck();
+  if (!hasFQ_) GaussCheck();
 
   infoStream_ << "========== Static solver " << std::endl;
   infoStream_ << *K_0_ << std::endl;