diff --git a/src/FFSetup.cpp b/src/FFSetup.cpp index fded99d3c..9f4388b08 100644 --- a/src/FFSetup.cpp +++ b/src/FFSetup.cpp @@ -19,6 +19,7 @@ along with this program, also can be found at // For Exotic style parameter header error checking #include +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 ¶m, 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 ¶m, 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 ¶m, 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 ¶m, 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 ¶m, 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 ¶m, 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 ¶m, 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 ¶m, 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 " diff --git a/src/FFSetup.h b/src/FFSetup.h index e7bd1f5a4..ca34608ed 100644 --- a/src/FFSetup.h +++ b/src/FFSetup.h @@ -40,7 +40,7 @@ class FFBase : public SearchableBase { std::string ReadKind(Reader ¶m, std::string const &firstKindName); std::vector &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 ¶m, 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 { diff --git a/src/Reader.h b/src/Reader.h index 28c4e28b9..5775114bf 100644 --- a/src/Reader.h +++ b/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;