- potentially fixed neighbor color stacking

[alexxy/gromacs-colorvec.git] / src / colorvec.cpp

diff --git a/src/colorvec.cpp b/src/colorvec.cpp
index 3fc5a3c9c056e9b1b698d7d993d00bb86b85218f..4acc2657bb4707a75b5fb512ccadb3b38774a9e0 100644 (file)
--- a/src/colorvec.cpp
+++ b/src/colorvec.cpp
@@ -289,13 +289,14 @@ void colorsStackingSearch(const t_pbc *inputPBC, const std::vector< gmx::RVec >
     gmx::RVec temp;
     #pragma omp parallel for ordered schedule(dynamic)
     for (size_t i = 0; i < colorsIndex.size(); ++i) {
-        for (size_t j {i + 1}; j < colorsIndex.size(); ++j) {
-            for (size_t k1 {0}; k1 < colorsIndex[i].size(); ++k1) {
-                for (size_t k2 {0}; k2 < colorsIndex[j].size(); ++k2) {
-                    pbc_dx(inputPBC, inputFrame[colorsIndex[i][k1].second], inputFrame[colorsIndex[j][k2].second], temp);
-                    if (temp.norm() <= radius) {
-                        ++outputTemp[i][j];
-                        ++outputTemp[j][i];
+        for (size_t j {0}; j < colorsIndex.size(); ++j) {
+            if (i != j) {
+                for (size_t k1 {0}; k1 < colorsIndex[i].size(); ++k1) {
+                    for (size_t k2 {0}; k2 < colorsIndex[j].size(); ++k2) {
+                        pbc_dx(inputPBC, inputFrame[colorsIndex[i][k1].second], inputFrame[colorsIndex[j][k2].second], temp);
+                        if (temp.norm() <= radius) {
+                            ++outputTemp[i][j];
+                        }
                     }
                 }
             }
@@ -308,10 +309,9 @@ void colorsStackingSearch(const t_pbc *inputPBC, const std::vector< gmx::RVec >
         output[i].resize(0);
     }
     for (size_t i {0}; i < outputTemp.size(); ++i) {
-        for (size_t j {i + 1}; j < outputTemp[i].size(); ++j) {
-            if (static_cast< float >(outputTemp[i][j]) >= static_cast< float >(colorsIndex.size()) * epsi) {
+        for (size_t j {0}; j < outputTemp[i].size(); ++j) {
+            if (static_cast< float >(outputTemp[i][j]) >= static_cast< float >(colorsIndex.size()) * epsi && i != j) {
                 output[i].push_back(j);
-                output[j].push_back(i);
             }
         }
     }
@@ -352,7 +352,7 @@ void anglesFileDump(const int frameNum, const std::string       &output,
     std::ofstream   file(output, std::ofstream::app);
     file << "frame =" << std::setw(8) << frameNum << std::endl;
     for (size_t i {0}; i < colorFormation.size(); ++i) {
-        file << "color #" << std::setw(3) << i;
+        file << "color #" << std::setw(3) << i + 1;
         if (toPeptide[i]) {
             file << std::setw(4) << "yes";
         } else {
@@ -387,7 +387,7 @@ void anglesFileDump(const int frameNum, const std::string       &output,
         if (inputStack[i].size() == 0) {
             file << std::setw(4) << "no []";
         } else {
-            file << std::setw(4) << "yes [ ";
+            file << std::setw(4) << "yes [";
             file << std::setw(4) << inputStack[i].size();
             for (size_t j {0}; j < inputStack[i].size(); ++j) {
                 file << std::setw(4) << inputStack[i][j];
@@ -401,12 +401,12 @@ void anglesFileDump(const int frameNum, const std::string       &output,
             for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
                 tempDBangle += dipoleBetaAngles[i][j];
             }
-            file << std::setw(8) << std::setprecision(2) << tempDBangle / dipoleBetaAngles[i].size();
+            file << std::setw(8) << std::setprecision(2) << tempDBangle / dipoleBetaAngles[i].size() << " [";
             for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
                 file << std::setw(8) << std::setprecision(2) << dipoleBetaAngles[i][j];
             }
         }
-        file << std::endl;
+        file << " ]" << std::endl;
     }
     file.close();
 }