Add data structures for FQ

Author: robertodr

src/interface/Input.cpp

@@ -166,12 +166,12 @@ void Input::reader(const std::string & filename) {
       std::vector<double> mono = chgdist.getDblVec("MONOPOLES");
       int j = 0;
       int n = int(mono.size() / 4);
-      multipoles_.monopoles = Eigen::VectorXd::Zero(n);
-      multipoles_.monopolesSites = Eigen::Matrix3Xd::Zero(3, n);
+      fragments_.monopoles = Eigen::VectorXd::Zero(n);
+      fragments_.monopolesSites = Eigen::Matrix3Xd::Zero(3, n);
       for (int i = 0; i < n; ++i) {
-        multipoles_.monopolesSites.col(i) =
+        fragments_.monopolesSites.col(i) =
             (Eigen::Vector3d() << mono[j], mono[j + 1], mono[j + 2]).finished();
-        multipoles_.monopoles(i) = mono[j + 3];
+        fragments_.monopoles(i) = mono[j + 3];
         j += 4;
       }
     }
@@ -180,16 +180,34 @@ void Input::reader(const std::string & filename) {
       std::vector<double> dipo = chgdist.getDblVec("DIPOLES");
       int j = 0;
       int n = int(dipo.size() / 6);
-      multipoles_.dipoles = Eigen::Matrix3Xd::Zero(3, n);
-      multipoles_.dipolesSites = Eigen::Matrix3Xd::Zero(3, n);
+      fragments_.dipoles = Eigen::Matrix3Xd::Zero(3, n);
+      fragments_.dipolesSites = Eigen::Matrix3Xd::Zero(3, n);
       for (int i = 0; i < n; ++i) {
-        multipoles_.dipolesSites.col(i) =
+        fragments_.dipolesSites.col(i) =
             (Eigen::Vector3d() << dipo[j], dipo[j + 1], dipo[j + 2]).finished();
-        multipoles_.dipoles.col(i) =
+        fragments_.dipoles.col(i) =
             (Eigen::Vector3d() << dipo[j + 3], dipo[j + 4], dipo[j + 5]).finished();
         j += 6;
       }
     }
+    // Set fluctuating charges
+    isFQ_ = false;
+    if (chgdist.getKey<std::vector<double> >("FQ").isDefined()) {
+      isFQ_ = true;
+      std::vector<double> fq = chgdist.getDblVec("FQ");
+      int j = 0;
+      int n = int(fq.size() / 5);
+      fragments_.FQSites = Eigen::Matrix3Xd::Zero(3, n);
+      fragments_.FQChi = Eigen::VectorXd::Zero(n);
+      fragments_.FQEta = Eigen::VectorXd::Zero(n);
+      for (int i = 0; i < n; ++i) {
+        fragments_.FQSites.col(i) =
+            (Eigen::Vector3d() << fq[j], fq[j + 1], fq[j + 2]).finished();
+        fragments_.FQChi(i) = fq[j + 3];
+        fragments_.FQEta(i) = fq[j + 4];
+        j += 5;
+      }
+    }
   }
 
   providedBy_ = std::string("API-side");
@@ -264,6 +282,7 @@ void Input::reader(const PCMInput & host_input) {
   correction_ = host_input.correction;
   hermitivitize_ = true;
   isDynamic_ = false;
+  isFQ_ = false;
 
   providedBy_ = std::string("host-side");
 }

src/interface/Input.hpp

@@ -116,6 +116,7 @@ class PCMSolver_EXPORT Input {
   bool hermitivitize() const { return hermitivitize_; }
   bool isDynamic() const { return isDynamic_; }
   double integratorScaling() const { return integratorScaling_; }
+  bool isFQ() const { return isFQ_; }
   /// @}
 
   /// Keeps track of who did the parsing: the API or the host program
@@ -130,7 +131,7 @@ class PCMSolver_EXPORT Input {
   BIOperatorData integratorParams() const;
   /// @}
 
-  ChargeDistribution multipoles() const { return multipoles_; }
+  ChargeDistribution fragments() const { return fragments_; }
   bool MEPfromMolecule() { return MEPfromMolecule_; }
 
   /// Operators
@@ -241,10 +242,12 @@ class PCMSolver_EXPORT Input {
   std::vector<double> origin_;
   /// Molecular geometry
   std::vector<double> geometry_;
+  /// Whether this is a FQ calculation
+  bool isFQ_;
   /// Whether to calculate the MEP from the molecular geometry
   bool MEPfromMolecule_;
-  /// Classical charge distribution of point multipoles
-  ChargeDistribution multipoles_;
+  /// Classical charge distribution
+  ChargeDistribution fragments_;
   /// Who performed the syntactic input parsing
   std::string providedBy_;
 };

src/utils/ChargeDistribution.hpp

@@ -52,6 +52,12 @@ struct ChargeDistribution {
   Eigen::Matrix3Xd dipoles;
   /*! Dipoles sites */
   Eigen::Matrix3Xd dipolesSites;
+  /*! FQ electronegativities */
+  Eigen::VectorXd FQChi;
+  /*! FQ hardnesses */
+  Eigen::VectorXd FQEta;
+  /*! FQ sites */
+  Eigen::Matrix3Xd FQSites;
 };
 
 /*! \typedef GFValue

tests/input/fq.inp

@@ -12,10 +12,10 @@ MEDIUM
   SOLVENT = WATER
 }
 
-FQ
+CHARGEDISTRIBUTION
 {
   # X, Y, Z, CHI, ETA
-  FRAGMENTS = [0.0, 0.0, 0.0, 1.0, 2.0,
-               0.0, 0.0, 0.0, 1.0, 2.0]
+  FQ = [0.0, 0.0, 0.0, 1.0, 2.0,
+        0.0, 0.0, 0.0, 1.0, 2.0]
 }
 

tools/pcmparser.py

@@ -361,12 +361,12 @@ def setup_keywords():
     # FQ section
     # Set a classical fluctuating charge force field
     # No additional spheres will be generated.
-    fq = getkw.Section('FQ', callback = verify_fq)
+    fq = Section('FQ', callback = verify_fq)
     # FQ force field
     # Valid values: array of doubles in format [x, y, z, chi, eta]
-    fq.add_kw('FRAGMENTS', 'DBL_ARRAY')
-    # TODO: smearing parameters
-    top.add_sect(fq)
+    charge_distribution.add_kw('FQ', 'DBL_ARRAY')
+    # TODO: smearing parameters for FQ
+    top.add_sect(charge_distribution)
 
     return top