Merge pull request #530 from GOMC-WSU/duplicate-dihedral

Detect and remove duplicate dihedrals with the same periodicity

src/FFSetup.cpp

@@ -19,6 +19,7 @@ along with this program, also can be found at
 // For Exotic style parameter header error checking
 #include <regex>
 
+const double EPSILON = 1.0e-4;
 const uint FFSetup::CHARMM_ALIAS_IDX = 0;
 const uint FFSetup::EXOTIC_ALIAS_IDX = 1;
 const std::string FFSetup::paramFileAlias[] = {"CHARMM-Style parameter file",
@@ -185,6 +186,12 @@ std::string FFBase::ReadKind(Reader &param, std::string const &firstKindName) {
     param.file >> tmp;
     merged += tmp;
   }
+  // Skip duplicates unless we allow multiple entries, such as for dihedrals
+  if (!multi && std::find(name.begin(), name.end(), merged) != name.end()) {
+    std::cout << "Warning: Ignoring duplicate entry for " << merged << " in "
+              << param.getFileName().c_str() << ". Using first entry.\n";
+  }
+  // Might insert duplicates but they will be ignored during execution
   if (!multi || std::find(name.begin(), name.end(), merged) == name.end())
     name.push_back(merged);
   return merged;
@@ -236,7 +243,7 @@ void Particle::Read(Reader &param, std::string const &firstVar) {
   // computation. But need the exponents to be large enough that the arithmetic
   // or geometric mean of any pair > 6.0. See FFParticle::Blend() for underlying
   // math.
-  double smallVal = 1e-20;
+  double smallVal = 1.0e-20;
   if (std::fabs(e) < smallVal) {
     e = 0.0;
     expN = 12.0; // Set to default (LJ) exponent.
@@ -246,12 +253,11 @@ void Particle::Read(Reader &param, std::string const &firstVar) {
     expN_1_4 = 12.0; // Set to default (LJ) exponent.
   }
   if ((expN - 6.0) < smallVal) {
-    std::cout << "ERROR: Mie exponent must be > 6!" << std::endl;
+    std::cout << "ERROR: Mie exponent must be > 6!\n";
     exit(EXIT_FAILURE);
   }
   if ((expN_1_4 - 6.0) < smallVal) {
-    std::cout << "ERROR: Mie exponent for 1-4 interactions must be > 6!"
-              << std::endl;
+    std::cout << "ERROR: Mie exponent for 1-4 interactions must be > 6!\n";
     exit(EXIT_FAILURE);
   }
 
@@ -380,14 +386,39 @@ void Dihedral::Read(Reader &param, std::string const &firstVar) {
   if (index == 0) {
     // set phase shift for n=0 to 90 degree
     // We will have C0 = Kchi (1 + cos(0 * phi + 90)) = Kchi
-    //this avoids double counting the C0 (constant offset) term, which is used force fields like TraPPE
+    // this avoids double counting the C0 (constant offset) term, which is used
+    // in force fields like TraPPE
     def = 90.00;
   }
-  Add(merged, coeff, index, def);
+  Add(param.getFileName(), merged, coeff, index, def);
   last = merged;
 }
-void Dihedral::Add(std::string const &merged, const double coeff,
-                   const uint index, const double def) {
+void Dihedral::Add(const std::string &fileName, const std::string &merged,
+                   const double coeff, const uint index, const double def) {
+  // Check for (and skip) duplicate periodicities for the same dihedral
+  // Generate an error and terminate if the duplicate dihedrals have different
+  // parameters
+  auto Kchi_it = Kchi[merged].begin();
+  auto delta_it = delta[merged].begin();
+  for (auto it = n[merged].begin(); it != n[merged].end(); ++it) {
+    // Found a duplicate dihedral
+    if (*it == index) {
+      if (std::fabs(*Kchi_it - EnConvIfCHARMM(coeff)) > EPSILON ||
+          std::fabs(*delta_it - geom::DegToRad(def)) > EPSILON) {
+        std::cout << "Error: Inconsistent dihedral parameters were found in "
+                  << fileName.c_str() << " for dihedral " << merged
+                  << " with periodicity " << index << "!\n";
+        exit(EXIT_FAILURE);
+      } else {
+        std::cout << "Warning: Ignoring duplicate periodicity of " << index
+                  << " for dihedral " << merged << " in " << fileName.c_str()
+                  << ".\n";
+        return;
+      }
+    }
+    Kchi_it++;
+    delta_it++;
+  }
   ++countTerms;
   Kchi[merged].push_back(EnConvIfCHARMM(coeff));
   n[merged].push_back(index);
@@ -399,7 +430,7 @@ void Improper::Read(Reader &param, std::string const &firstVar) {
   uint index;
   std::string merged = ReadKind(param, firstVar);
   // If new value
-  if (validname(merged) == true) {
+  if (validname(merged)) {
     param.file >> coeff >> index >> def;
     if (!param.file.good()) {
       std::cout << "Error: Incomplete Improper parameters was found in "
@@ -420,7 +451,7 @@ void CMap::Read(Reader &param, std::string const &firstVar) {
   uint index;
   std::string merged = ReadKind(param, firstVar);
   // If new value
-  if (validname(merged) == true) {
+  if (validname(merged)) {
     param.file >> coeff >> index >> def;
     if (!param.file.good()) {
       std::cout << "Error: Incomplete Improper parameters was found in "
@@ -430,19 +461,19 @@ void CMap::Read(Reader &param, std::string const &firstVar) {
     Add(coeff, def);
   }
 }
-// Currently dummy method, exact same as improper
+// Currently dummy method, exactly the same as improper
 void CMap::Add(const double coeff, const double def) {
   Komega.push_back(EnConvIfCHARMM(coeff));
   omega0.push_back(def);
 }
 
-// Currently dummy method, exact same as improper
+// Currently dummy method, exactly the same as improper
 void HBond::Read(Reader &param, std::string const &firstVar) {
   double coeff, def;
   uint index;
   std::string merged = ReadKind(param, firstVar);
   // If new value
-  if (validname(merged) == true) {
+  if (validname(merged)) {
     param.file >> coeff >> index >> def;
     if (!param.file.good()) {
       std::cout << "Error: Incomplete Improper parameters was found in "

src/FFSetup.h

@@ -40,7 +40,7 @@ class FFBase : public SearchableBase {
   std::string ReadKind(Reader &param, std::string const &firstKindName);
   std::vector<std::string> &getnamelist() { return name; }
 
-  bool validname(std::string &merged) {
+  bool validname(const std::string &merged) const {
     return (std::count(name.begin(), name.end(), merged) > 0);
   }
   std::string getname(const uint i) const { return name[i]; }
@@ -146,8 +146,8 @@ class Dihedral : public ReadableBaseWithFirst, public FFBase {
 public:
   Dihedral(void) : FFBase(4, true), last(""), countTerms(0) {}
   virtual void Read(Reader &param, std::string const &firstVar);
-  void Add(std::string const &merged, const double coeff, const uint index,
-           const double def);
+  void Add(const std::string &fileName, const std::string &merged,
+           const double coeff, const uint index, const double def);
   uint getTerms() const { return countTerms; }
   uint append(std::string &s, double *Kchi_in, double *delta_in, uint *n_in,
               uint count) const {

src/Reader.h

@@ -32,6 +32,7 @@ inline std::ifstream &operator>>(std::ifstream &is, bool &val) {
 
 class Reader {
 public:
+  friend inline std::ifstream &operator>>(std::ifstream &is, bool &val);
   Reader(std::string const &name, std::string const &alias,
          bool useSkipW = false, std::string const *const skipW = NULL,
          bool useSkipC = false, std::string const *const skipC = NULL,
@@ -98,6 +99,7 @@ class Reader {
   }
 
   bool Read(std::string &firstItem);
+  std::string getFileName() const { return fileName; };
 
   // Make public to expose object.
   std::ifstream file;