Make it possible to run nonpolarizable MM

Note this is a horrible, horrible, horrible HACK

Author: robertodr

src/interface/Input.cpp

@@ -198,6 +198,7 @@ void Input::reader(const std::string & filename) {
   const Section & mmfq = input_.getSect("MMFQ");
   if (mmfq.isDefined()) {
     isFQ_ = true;
+    isNonPolarizable_ = mmfq.getBool("NONPOLARIZABLE");
     fragments_.nSitesPerFragment =
         static_cast<PCMSolverIndex>(mmfq.getInt("SITESPERFRAGMENT"));
     std::vector<double> sites = mmfq.getDblVec("SITES");

src/interface/Input.hpp

@@ -119,6 +119,7 @@ class PCMSolver_EXPORT Input {
   bool isDynamic() const { return isDynamic_; }
   double integratorScaling() const { return integratorScaling_; }
   bool isFQ() const { return isFQ_; }
+  bool isNonPolarizable() const { return isNonPolarizable_; }
   /// @}
 
   /// Keeps track of who did the parsing: the API or the host program
@@ -248,6 +249,8 @@ class PCMSolver_EXPORT Input {
   std::vector<double> geometry_;
   /// Whether this is a FQ calculation
   bool isFQ_;
+  /// Whether this is a nonpolarizable MM calculation
+  bool isNonPolarizable_;
   /// Whether to calculate the MEP from the molecular geometry
   bool MEPfromMolecule_;
   /// Whether to calculate the MEP from the charge distribution

src/interface/Meddle.cpp

@@ -155,16 +155,16 @@ void Meddle::CTORBody() {
     TIMER_ON("Meddle::initCavity");
     initCavity();
     TIMER_OFF("Meddle::initCavity");
+  }
 
-    TIMER_ON("Meddle::initStaticSolver");
-    initStaticSolver();
-    TIMER_OFF("Meddle::initStaticSolver");
+  TIMER_ON("Meddle::initStaticSolver");
+  initStaticSolver();
+  TIMER_OFF("Meddle::initStaticSolver");
 
-    if (input_.isDynamic()) {
-      TIMER_ON("Meddle::initDynamicSolver");
-      initDynamicSolver();
-      TIMER_OFF("Meddle::initDynamicSolver");
-    }
+  if (input_.isDynamic()) {
+    TIMER_ON("Meddle::initDynamicSolver");
+    initDynamicSolver();
+    TIMER_OFF("Meddle::initDynamicSolver");
   }
 }
 
@@ -335,10 +335,12 @@ double pcm::Meddle::computePolarizationEnergy(const std::string & mep_name,
     // Dot product of MEP + Electronegativities times Fluctuating charges
     energy = (functions_.at(mep_name) + input_.fragments().chi)
                  .dot(functions_.at(asc_name));
+    if (input_.isNonPolarizable()) { /* HACK Nonpolarizable doesn't need 1/2 */
+      energy *= 2.0;
+    }
   } else { /* Pure PCM calculation */
     // Dot product of MEP and ASC surface function
-    energy =
-        functions_.at(mep_name).dot(functions_.at(asc_name));
+    energy = functions_.at(mep_name).dot(functions_.at(asc_name));
   }
   return (energy / 2.0);
 }
@@ -372,8 +374,13 @@ void pcm::Meddle::computeASC(const std::string & mep_name,
   // Get the proper iterators
   SurfaceFunctionMapConstIter iter_pot = functions_.find(mep_name);
   Eigen::VectorXd asc = Eigen::VectorXd::Zero(iter_pot->second.size());
-  if (hasFQ_) { /* MMFQ calculation */
-    asc = FQ_->computeCharge(iter_pot->second);
+  if (hasFQ_) {                      /* MMFQ calculation */
+    if (input_.isNonPolarizable()) { /* HACK We store point charges in the eta vector
+                                      */
+      asc = input_.fragments().eta;
+    } else {
+      asc = FQ_->computeCharge(iter_pot->second);
+    }
   } else { /* Pure PCM calculation */
     asc = K_0_->computeCharge(iter_pot->second, irrep);
     // Renormalize
@@ -402,9 +409,11 @@ void pcm::Meddle::computeResponseASC(const std::string & mep_name,
   // Get the proper iterators
   SurfaceFunctionMapConstIter iter_pot = functions_.find(mep_name);
   Eigen::VectorXd asc = Eigen::VectorXd::Zero(iter_pot->second.size());
-  if (hasFQ_) { /* MMFQ calculation */
-    // Do NOT add classical (electronegativities) contributions to RHS
-    asc = FQ_->computeCharge(iter_pot->second, false);
+  if (hasFQ_) {                       /* MMFQ calculation */
+    if (!input_.isNonPolarizable()) { /* HACK Can we do it more cleanly/clearly? */
+      // Do NOT add classical (electronegativities) contributions to RHS
+      asc = FQ_->computeCharge(iter_pot->second, false);
+    }
   } else { /* Pure PCM calculation */
     if (hasDynamic_) {
       asc = K_d_->computeCharge(iter_pot->second, irrep);
@@ -726,7 +735,8 @@ void Meddle::initStaticSolver() {
   delete biop;
 
   // Perform Gauss' theorem check for nuclear charges
-  if (!hasFQ_) GaussCheck();
+  if (!hasFQ_)
+    GaussCheck();
 
   infoStream_ << "========== Static solver " << std::endl;
   infoStream_ << *K_0_ << std::endl;
@@ -759,8 +769,8 @@ void Meddle::initDynamicSolver() {
 }
 
 void Meddle::initMMFQ() {
-  FQ_ = new mmfq::FQOhno(input_.fragments());
-  size_ = pcm::make_tuple(input_.fragments().sites.cols(),
+  FQ_ = new mmfq::FQOhno(input_.fragments(), input_.isNonPolarizable());
+  size_ = std::make_tuple(input_.fragments().sites.cols(),
                           input_.fragments().sites.cols());
   hasFQ_ = true;
 

src/mmfq/CMakeLists.txt

@@ -1,33 +1,12 @@
+target_sources(pcm-objlib
+  PRIVATE
+    ${CMAKE_CURRENT_SOURCE_DIR}/FQOhno.cpp
+  )
+
 # List of headers
 list(APPEND headers_list
   FQOhno.hpp
   )
-
-# List of sources
-list(APPEND sources_list
-  FQOhno.cpp
-  )
-
-add_library(mmfq OBJECT ${sources_list} ${headers_list})
-set_target_properties(mmfq
-  PROPERTIES
-    POSITION_INDEPENDENT_CODE 1
-    CXX_VISIBILITY_PRESET hidden
-    VISIBILITY_INLINES_HIDDEN 1
-  )
-target_compile_definitions(mmfq
-  PUBLIC
-    PCMSolver_EXPORTS
-  )
-if(BUILD_CUSTOM_BOOST)
-  add_dependencies(mmfq custom_boost)
-endif()
-add_dependencies(mmfq generate-config-hpp)
-target_compile_options(mmfq
-  PRIVATE
-    "$<$<CONFIG:DEBUG>:${EXDIAG_CXX_FLAGS}>"
-  )
-
 # Sets install directory for all the headers in the list
 foreach(_header ${headers_list})
     install(FILES ${_header} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}/mmfq)

src/mmfq/FQOhno.cpp

@@ -34,10 +34,6 @@
 
 namespace pcm {
 namespace mmfq {
-using utils::MMFQ;
-
-FQOhno::FQOhno(const MMFQ & ff) : mmfq_(ff) { buildSystemMatrix_impl(); }
-
 void FQOhno::buildSystemMatrix_impl() {
   PCMSolverIndex nFragments = mmfq_.nFragments;
   PCMSolverIndex nSitesPerFragment = mmfq_.nSitesPerFragment;
@@ -80,6 +76,7 @@ Eigen::VectorXd FQOhno::computeCharge_impl(const Eigen::VectorXd & potential,
 
 std::ostream & FQOhno::printSolver(std::ostream & os) {
   os << "Fluctuating charge solver type: Ohno" << std::endl;
+  if (nonPolarizable_) os << "Nonpolarizable force field" << std::endl;
   os << "Number of fragments = " << mmfq_.nFragments << std::endl;
   os << "Number of sites per fragment = " << mmfq_.nSitesPerFragment << std::endl;
   os << "Number of sites = " << mmfq_.nFragments * mmfq_.nSitesPerFragment;

src/mmfq/FQOhno.hpp

@@ -43,15 +43,19 @@ namespace mmfq {
  *  to solve for the charges.
  *  This avoids computing and storing the inverse explicitly.
  */
-class FQOhno __final {
+class FQOhno final {
 public:
-  FQOhno() {}
+  FQOhno() = default;
   /*! \brief Construct solver
    *  \param[in] ff the fluctuating charges force field
+   *  \param[in] nonpol whether this is a nonpolarizable MM calculation
    */
-  FQOhno(const utils::MMFQ & ff);
+  FQOhno(const utils::MMFQ & ff, bool nonpol = false)
+      : nonPolarizable_(nonpol), built_(false), mmfq_(ff) {
+    if (!nonPolarizable_)
+      buildSystemMatrix_impl();
+  }
 
-  ~FQOhno() {}
   friend std::ostream & operator<<(std::ostream & os, FQOhno & solver) {
     return solver.printSolver(os);
   }
@@ -63,6 +67,7 @@ class FQOhno __final {
   }
 
 private:
+  bool nonPolarizable_;
   bool built_;
   utils::MMFQ mmfq_;
   /*! D_\lambda matrix, not symmetry blocked */

src/utils/CMakeLists.txt

@@ -18,6 +18,7 @@ list(APPEND headers_list
   Logger.hpp
   LoggerImpl.hpp
   MathUtils.hpp
+  MMFQ.hpp
   Molecule.hpp
   QuadratureRules.hpp
   Solvent.hpp

tests/mmfq/CMakeLists.txt

@@ -1,12 +1,14 @@
-add_library(mmfq-tests OBJECT mmfq-ohno.cpp)
-if(BUILD_CUSTOM_BOOST)
-  add_dependencies(mmfq-tests custom_boost)
-endif()
-
-# Copy reference files to ${PROJECT_BINARY_DIR}/tests/mmfq (aka ${CMAKE_CURRENT_BINARY_DIR})
-list(APPEND reference_files FQ-2xH2O.npy)
-file(COPY ${reference_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
-
-# mmfq.cpp test
-add_Catch_test(mmfq-ohno "mmfq;ohno")
+target_sources(unit_tests
+  PRIVATE
+    ${CMAKE_CURRENT_SOURCE_DIR}/mmfq-ohno.cpp
+  )
 
+add_Catch_test(
+  NAME
+    mmfq-ohno
+  LABELS
+    mmfq
+    ohno
+  REFERENCE_FILES
+    FQ-2xH2O.npy
+  )

tools/pcmparser.py

@@ -362,6 +362,9 @@ def setup_keywords():
     # Set a classical fluctuating charge force field
     # No additional spheres will be generated.
     mmfq = Section('MMFQ', callback = verify_mmfq)
+    # Valid values: boolean
+    # Default: False
+    mmfq.add_kw('NONPOLARIZABLE', 'BOOL', False)
     # Valid values: integer
     mmfq.add_kw('SITESPERFRAGMENT', 'INT', 3)
     # FQ force field
@@ -583,6 +586,7 @@ def verify_mmfq(section):
         print('Empty or incoherent SITES array')
         sys.exit(1)
 
+
 def verify_spheres(keyword):
     length = len(keyword.get())
     if (length % 4 != 0):