Avoid adding classical portion to RHS when doing response

Author: robertodr

src/interface/Meddle.cpp

@@ -402,7 +402,8 @@ void pcm::Meddle::computeResponseASC(const std::string & mep_name,
   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);
+    // 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);

src/mmfq/FQOhno.cpp

@@ -68,10 +68,13 @@ void FQOhno::buildSystemMatrix_impl() {
   built_ = true;
 }
 
-Eigen::VectorXd FQOhno::computeCharge_impl(const Eigen::VectorXd & potential) const {
-  // We have to add electronegativities on top of the potential
+Eigen::VectorXd FQOhno::computeCharge_impl(const Eigen::VectorXd & potential,
+                                           bool scf) const {
   Eigen::VectorXd RHS = Eigen::VectorXd::Zero(Dlambda_.rows());
-  RHS.head(potential.size()) = potential + mmfq_.chi;
+  RHS.head(potential.size()) = potential;
+  // If doing SCF, we have to add electronegativities on top of the potential
+  if (scf)
+    RHS.head(potential.size()) += mmfq_.chi;
   return -Dlambda_.ldlt().solve(RHS).head(potential.size());
 }
 

src/mmfq/FQOhno.hpp

@@ -57,10 +57,11 @@ class FQOhno __final {
   friend std::ostream & operator<<(std::ostream & os, FQOhno & solver) {
     return solver.printSolver(os);
   }
-  Eigen::VectorXd computeCharge(const Eigen::VectorXd & potential) const {
+  Eigen::VectorXd computeCharge(const Eigen::VectorXd & potential,
+                                bool scf = true) const {
     if (!built_)
       PCMSOLVER_ERROR("MMFQ matrix not calculated yet");
-    return computeCharge_impl(potential);
+    return computeCharge_impl(potential, scf);
   }
 
 private:
@@ -70,7 +71,8 @@ class FQOhno __final {
   Eigen::MatrixXd Dlambda_;
 
   void buildSystemMatrix_impl();
-  Eigen::VectorXd computeCharge_impl(const Eigen::VectorXd & potential) const;
+  Eigen::VectorXd computeCharge_impl(const Eigen::VectorXd & potential,
+                                     bool scf = true) const;
   std::ostream & printSolver(std::ostream & os);
 };
 } // namespace mmfq