From f42c704fde2f66f7bd44096e68bdcfd39cc1ab88 Mon Sep 17 00:00:00 2001 From: msoroush Date: Thu, 26 Sep 2019 16:35:48 -0400 Subject: [PATCH 1/4] Fix to the issue #142 --- src/MolSetup.cpp | 227 ++++++++++++++++++++++------------------------- 1 file changed, 106 insertions(+), 121 deletions(-) diff --git a/src/MolSetup.cpp b/src/MolSetup.cpp index 5148cb0e2..ab6d14c63 100644 --- a/src/MolSetup.cpp +++ b/src/MolSetup.cpp @@ -70,17 +70,20 @@ int ReadPSFAtoms(FILE* psf, //pre: stream is before !BONDS post: stream is in bond section just after //the first appearance of the last molecule int ReadPSFBonds(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom); + std::vector >& firstAtom, + const uint nbonds); //adds angles in psf to kindMap //pre: stream is before !NTHETA post: stream is in angle section just //after the first appearance of the last molecule int ReadPSFAngles(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom); + std::vector >& firstAtom, + const uint nangles); //adds dihedrals in psf to kindMap //pre: stream is before !NPHI post: stream is in dihedral section just //after the first appearance of the last molecule int ReadPSFDihedrals(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom); + std::vector >& firstAtom, + const uint ndihedrals); } @@ -567,65 +570,61 @@ int ReadPSF(const char* psfFilename, MolMap& kindMap) //build list of start particles for each type, so we can find it and skip //everything else std::vector > firstAtomLookup; - for (MolMap::iterator it = kindMap.begin(); - it != kindMap.end(); ++it) { - firstAtomLookup.push_back(std::make_pair(it->second.firstAtomID, it->first)); + for (MolMap::iterator it = kindMap.begin(); it != kindMap.end(); ++it) { + firstAtomLookup.push_back(std::make_pair(it->second.firstAtomID,it->first)); } std::sort(firstAtomLookup.begin(), firstAtomLookup.end()); //find bond header+count while (strstr(input, "!NBOND") == NULL) { check = fgets(input, 511, psf); if (check == NULL) { - fprintf(stderr, "ERROR: Unable to read bonds from PSF file %s", psfFilename); + fprintf(stderr, "ERROR: Unable to read bonds from PSF file %s", + psfFilename); fclose(psf); return READERROR; } } //make sure molecule has bonds, appears before !NBOND count = atoi(input); - if (count != 0) { - if (ReadPSFBonds(psf, kindMap, firstAtomLookup) == READERROR) { - fclose(psf); - return READERROR; - } + if (ReadPSFBonds(psf, kindMap, firstAtomLookup, count) == READERROR) { + fclose(psf); + return READERROR; } //find angle header+count fseek(psf, 0, SEEK_SET); while (strstr(input, "!NTHETA") == NULL) { check = fgets(input, 511, psf); if (check == NULL) { - fprintf(stderr, "ERROR: Unable to read angles from PSF file %s", psfFilename); + fprintf(stderr, "ERROR: Unable to read angles from PSF file %s", + psfFilename); fclose(psf); return READERROR; } } //make sure molecule has angles, count appears before !NTHETA count = atoi(input); - if (count != 0) { - if (ReadPSFAngles(psf, kindMap, firstAtomLookup) == READERROR) { - fclose(psf); - return READERROR; - } + if (ReadPSFAngles(psf, kindMap, firstAtomLookup, count) == READERROR) { + fclose(psf); + return READERROR; } //find dihedrals header+count fseek(psf, 0, SEEK_SET); while (strstr(input, "!NPHI") == NULL) { check = fgets(input, 511, psf); if (check == NULL) { - fprintf(stderr, "ERROR: Unable to read dihedrals from PSF file %s", psfFilename); + fprintf(stderr, "ERROR: Unable to read dihedrals from PSF file %s", + psfFilename); fclose(psf); return READERROR; } } //make sure molecule has dihs, count appears before !NPHI count = atoi(input); - if (count != 0) { - if (ReadPSFDihedrals(psf, kindMap, firstAtomLookup) == READERROR) { - fclose(psf); - return READERROR; - } + if (ReadPSFDihedrals(psf, kindMap, firstAtomLookup, count) == READERROR) { + fclose(psf); + return READERROR; } - + fclose(psf); return nAtoms; @@ -682,33 +681,34 @@ int ReadPSFAtoms(FILE* psf, MolMap& kindMap, unsigned int nAtoms) //pre: stream is before !BONDS post: stream is in bond section just after //the first appearance of the last molecule int ReadPSFBonds(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom) + std::vector >& firstAtom, + const uint nbonds) { unsigned int atom0, atom1; - int dummy = fscanf(psf, "%u %u", &atom0, &atom1); - UNUSED(dummy); - for (unsigned int i = 0; i < firstAtom.size(); ++i) { - MolKind& currentMol = kindMap[firstAtom[i].second]; - //continue if atom has no bonds - if (currentMol.atoms.size() < 2) - continue; - //index of first atom in moleule - unsigned int molBegin = firstAtom[i].first; - //index AFTER last atom in molecule - unsigned int molEnd = molBegin + currentMol.atoms.size(); - while (atom0 < molBegin || atom0 >= molEnd) { - dummy = fscanf(psf, "%u %u", &atom0, &atom1); - if (feof(psf) || ferror(psf)) { - fprintf(stderr, "ERROR: Could not find all bonds in PSF file "); - return READERROR; - } + int dummy; + for (uint n = 0; n < nbonds; n++) { + dummy = fscanf(psf, "%u %u", &atom0, &atom1); + if(dummy != 2) { + fprintf(stderr, "ERROR: Incorrect Number of bonds in PSF file "); + return READERROR; + } else if (feof(psf) || ferror(psf)) { + fprintf(stderr, "ERROR: Could not find all bonds in PSF file "); + return READERROR; } - //read in bonds - while (atom0 >= molBegin && atom0 < molEnd) { - currentMol.bonds.push_back(Bond(atom0 - molBegin, atom1 - molBegin)); - dummy = fscanf(psf, "%u %u", &atom0, &atom1); - if(dummy != 2) - break; + + //loop to find the molecule kind with this bond + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + //index of first atom in molecule + unsigned int molBegin = firstAtom[i].first; + //index AFTER last atom in molecule + unsigned int molEnd = molBegin + currentMol.atoms.size(); + //assign the bond + if (atom0 >= molBegin && atom0 < molEnd) { + currentMol.bonds.push_back(Bond(atom0 - molBegin, atom1 - molBegin)); + //once we found the molecule kind, break from the loop + break; + } } } return 0; @@ -718,34 +718,35 @@ int ReadPSFBonds(FILE* psf, MolMap& kindMap, //pre: stream is before !NTHETA post: stream is in angle section just after //the first appearance of the last molecule int ReadPSFAngles(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom) + std::vector >& firstAtom, + const uint nangles) { unsigned int atom0, atom1, atom2; - int dummy = fscanf(psf, "%u %u %u", &atom0, &atom1, &atom2); - UNUSED(dummy); - for (unsigned int i = 0; i < firstAtom.size(); ++i) { - MolKind& currentMol = kindMap[firstAtom[i].second]; - //continue if atom has no angles - if (currentMol.atoms.size() < 3) - continue; - //index of first atom in moleule - unsigned int molBegin = firstAtom[i].first; - //index AFTER last atom in molecule - unsigned int molEnd = molBegin + currentMol.atoms.size(); - while (atom0 < molBegin || atom0 >= molEnd) { - dummy = fscanf(psf, "%u %u %u", &atom0, &atom1, &atom2); - if (feof(psf) || ferror(psf)) { - fprintf(stderr, "ERROR: Could not find all angles in PSF file "); - return READERROR; - } + int dummy; + for (uint n = 0; n < nangles; n++) { + dummy = fscanf(psf, "%u %u %u", &atom0, &atom1, &atom2); + if(dummy != 3) { + fprintf(stderr, "ERROR: Incorrect Number of angles in PSF file "); + return READERROR; + } else if (feof(psf) || ferror(psf)) { + fprintf(stderr, "ERROR: Could not find all angles in PSF file "); + return READERROR; } - //read in angles - while (atom0 >= molBegin && atom0 < molEnd) { - currentMol.angles.push_back(Angle(atom0 - molBegin, atom1 - molBegin, - atom2 - molBegin)); - dummy = fscanf(psf, "%u %u %u", &atom0, &atom1, &atom2); - if(dummy != 3) - break; + + //loop to find the molecule kind with this angles + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + //index of first atom in moleule + unsigned int molBegin = firstAtom[i].first; + //index AFTER last atom in molecule + unsigned int molEnd = molBegin + currentMol.atoms.size(); + //assign the angle + if (atom0 >= molBegin && atom0 < molEnd) { + currentMol.angles.push_back(Angle(atom0 - molBegin, atom1 - molBegin, + atom2 - molBegin)); + //once we found the molecule kind, break from the loop + break; + } } } return 0; @@ -772,60 +773,44 @@ bool ContainsDihedral(const std::vector& vec, const int dih[]) //the first appearance of the last molecule // int ReadPSFDihedrals(FILE* psf, MolMap& kindMap, - std::vector >& firstAtom) + std::vector >& firstAtom, const uint ndihedrals) { Dihedral dih(0, 0, 0, 0); - - int dummy = fscanf(psf, "%u %u %u %u", &dih.a0, &dih.a1, &dih.a2, &dih.a3); - UNUSED(dummy); - //for all atoms - for (unsigned int i = 0; i < firstAtom.size(); ++i) { - MolKind& currentMol = kindMap[firstAtom[i].second]; - //continue if molecule has no dihedrals - if (currentMol.atoms.size() < 4) - continue; - //index of first atom in moleule - unsigned int molBegin = firstAtom[i].first; - //index AFTER last atom in molecule - unsigned int molEnd = molBegin + currentMol.atoms.size(); - //continue if molecule has more that 3 atoms but has no dihedrals - if(i == 0) { - // if it is the first molecule and index of dihedral is greater than - // molBegin, it means it does not have any dihedral. It works when - // we have only two molecule kinds. - if(dih.a0 > molBegin && dih.a0 > molEnd) { - continue; - } + int dummy; + for (uint n = 0; n < ndihedrals; n++) { + dummy = fscanf(psf, "%u %u %u %u", &dih.a0, &dih.a1, &dih.a2, &dih.a3); + if(dummy != 4) { + fprintf(stderr, "ERROR: Incorrect Number of dihedrals in PSF file "); + return READERROR; + } else if (feof(psf) || ferror(psf)) { + fprintf(stderr, "ERROR: Could not find all dihedrals in PSF file "); + return READERROR; } - //scan to to first appearance of molecule - while (dih.a0 < molBegin || dih.a0 >= molEnd) { - dummy = fscanf(psf, "%u %u %u %u", &dih.a0, &dih.a1, &dih.a2, &dih.a3); - if (feof(psf) || ferror(psf)) { - fprintf(stderr, "ERROR: Could not find all dihedrals in PSF file "); - return READERROR; - } - //for case that molecule has more thatn 3 atoms and represend as - //second molecule but it has no dihedral. It works when - // we have only two molecule kinds and no improper - if (dih.a0 == 0) + + //loop to find the molecule kind with this dihedral + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + //index of first atom in moleule + unsigned int molBegin = firstAtom[i].first; + //index AFTER last atom in molecule + unsigned int molEnd = molBegin + currentMol.atoms.size(); + //assign dihedral + if (dih.a0 >= molBegin && dih.a0 < molEnd) { + dih.a0 -= molBegin; + dih.a1 -= molBegin; + dih.a2 -= molBegin; + dih.a3 -= molBegin; + //some xplor PSF files have duplicate dihedrals, we need to ignore these + if (std::find(currentMol.dihedrals.begin(), currentMol.dihedrals.end(), + dih) == currentMol.dihedrals.end()) { + currentMol.dihedrals.push_back(dih); + } + //once we found the molecule kind, break from the loop break; - } - //read in dihedrals - while (dih.a0 >= molBegin && dih.a0 < molEnd) { - dih.a0 -= molBegin; - dih.a1 -= molBegin; - dih.a2 -= molBegin; - dih.a3 -= molBegin; - //some xplor PSF files have duplicate dihedrals, we need to ignore these - if (std::find(currentMol.dihedrals.begin(), currentMol.dihedrals.end(), - dih) == currentMol.dihedrals.end()) { - currentMol.dihedrals.push_back(dih); } - dummy = fscanf(psf, "%u %u %u %u", &dih.a0, &dih.a1, &dih.a2, &dih.a3); - if(dummy != 4) - break; } } return 0; } + } From f9f6cc133fc79fc5153bd98cdb2bccfd2d90fbf1 Mon Sep 17 00:00:00 2001 From: msoroush Date: Thu, 26 Sep 2019 17:05:29 -0400 Subject: [PATCH 2/4] Generate warning if we expct bond, angle, or dihedral for each molecule kind. Fix to issue #142 --- src/MolSetup.cpp | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/MolSetup.cpp b/src/MolSetup.cpp index ab6d14c63..ac79b09b0 100644 --- a/src/MolSetup.cpp +++ b/src/MolSetup.cpp @@ -685,7 +685,8 @@ int ReadPSFBonds(FILE* psf, MolMap& kindMap, const uint nbonds) { unsigned int atom0, atom1; - int dummy; + int dummy; + std::vector defined(firstAtom.size(), false); for (uint n = 0; n < nbonds; n++) { dummy = fscanf(psf, "%u %u", &atom0, &atom1); if(dummy != 2) { @@ -707,10 +708,19 @@ int ReadPSFBonds(FILE* psf, MolMap& kindMap, if (atom0 >= molBegin && atom0 < molEnd) { currentMol.bonds.push_back(Bond(atom0 - molBegin, atom1 - molBegin)); //once we found the molecule kind, break from the loop + defined[i] = true; break; } } } + //Check if we defined all bonds + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + if(currentMol.atoms.size() > 1 && !defined[i]) { + std::cout << "Warning: Bond is missing for " << firstAtom[i].second + << " !\n"; + } + } return 0; } @@ -723,6 +733,7 @@ int ReadPSFAngles(FILE* psf, MolMap& kindMap, { unsigned int atom0, atom1, atom2; int dummy; + std::vector defined(firstAtom.size(), false); for (uint n = 0; n < nangles; n++) { dummy = fscanf(psf, "%u %u %u", &atom0, &atom1, &atom2); if(dummy != 3) { @@ -745,10 +756,19 @@ int ReadPSFAngles(FILE* psf, MolMap& kindMap, currentMol.angles.push_back(Angle(atom0 - molBegin, atom1 - molBegin, atom2 - molBegin)); //once we found the molecule kind, break from the loop + defined[i] = true; break; } } } + //Check if we defined all angles + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + if(currentMol.atoms.size() > 2 && !defined[i]) { + std::cout << "Warning: Angle is missing for " << firstAtom[i].second + << " !\n"; + } + } return 0; } @@ -777,6 +797,7 @@ int ReadPSFDihedrals(FILE* psf, MolMap& kindMap, { Dihedral dih(0, 0, 0, 0); int dummy; + std::vector defined(firstAtom.size(), false); for (uint n = 0; n < ndihedrals; n++) { dummy = fscanf(psf, "%u %u %u %u", &dih.a0, &dih.a1, &dih.a2, &dih.a3); if(dummy != 4) { @@ -806,10 +827,19 @@ int ReadPSFDihedrals(FILE* psf, MolMap& kindMap, currentMol.dihedrals.push_back(dih); } //once we found the molecule kind, break from the loop + defined[i] = true; break; } } } + //Check if we defined all dihedrals + for (unsigned int i = 0; i < firstAtom.size(); ++i) { + MolKind& currentMol = kindMap[firstAtom[i].second]; + if(currentMol.atoms.size() > 3 && !defined[i]) { + std::cout << "Warning: Dihedral is missing for " << firstAtom[i].second + << " !\n"; + } + } return 0; } From 0c40b766d224202b2aee5800eed55a77daeb9dd4 Mon Sep 17 00:00:00 2001 From: msoroush Date: Sun, 6 Oct 2019 13:58:38 -0400 Subject: [PATCH 3/4] Suggested fix to the issue #143 and #123 --- src/CellList.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CellList.h b/src/CellList.h index 88c9748b3..09da998e9 100644 --- a/src/CellList.h +++ b/src/CellList.h @@ -85,6 +85,12 @@ inline int CellList::PositionToCell(const XYZ& posRef, int box) const int x = (int)(pos.x / cellSize[box].x); int y = (int)(pos.y / cellSize[box].y); int z = (int)(pos.z / cellSize[box].z); + //Check the cell number to avoid segfult for coordinates close to axis + //x, y, and z should never be equal or greater than number of cells in x, y, + // and z axis, respectively. + x -= (x == edgeCells[box][0] ? 1 : 0); + y -= (y == edgeCells[box][1] ? 1 : 0); + z -= (z == edgeCells[box][2] ? 1 : 0); return x * edgeCells[box][1] * edgeCells[box][2] + y * edgeCells[box][2] + z; } From 7135ecb5e76d2cb444c424b357067bbda3d5cb86 Mon Sep 17 00:00:00 2001 From: msoroush Date: Sun, 6 Oct 2019 14:04:12 -0400 Subject: [PATCH 4/4] Fix printing config file info for MEMC parameters --- src/ConfigSetup.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ConfigSetup.cpp b/src/ConfigSetup.cpp index d0a218a8a..1a3887c49 100644 --- a/src/ConfigSetup.cpp +++ b/src/ConfigSetup.cpp @@ -339,7 +339,7 @@ void ConfigSetup::Init(const char *fileName) exit(EXIT_FAILURE); } if(line.size() >= 3) { - printf("%-41s", "Info: Atom Names in BackBone of Small Molecule Kind"); + printf("%-41s", "Info: Atom Names in BackBone of Small Molecule Kind "); for(uint i = 1; i < line.size() - 1; i += 2) { if(i != 1) { printf(" , "); @@ -362,7 +362,7 @@ void ConfigSetup::Init(const char *fileName) exit(EXIT_FAILURE); } if(line.size() >= 3) { - printf("%-41s", "Info: Atom Names in BackBone of Large Molecule Kind"); + printf("%-41s", "Info: Atom Names in BackBone of Large Molecule Kind "); for(uint i = 1; i < line.size() - 1; i += 2) { if(i != 1) { printf(" , ");