src/domains.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43 #include <iostream>
  44 #include <chrono>
  45 #include <cstdint>
  46 #include <cfloat>
  47 #include <omp.h>
  48
  49 #include <gromacs/trajectoryanalysis.h>
  50 #include <gromacs/pbcutil/pbc.h>
  51 #include <gromacs/utility/smalloc.h>
  52 #include <gromacs/math/do_fit.h>
  53
  54 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/topology/index.h"
  57
  58 //#include "/home/titov_ai/Develop/fittingLib/fittinglib.h"
  59
  60 #define MAX_NTRICVEC 12
  61
  62 using namespace gmx;
  63
  64 using gmx::RVec;
  65
  66 double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
  67     return  sqrt(   pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
  68                     pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
  69                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  70 }
  71
  72 void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
  73     double t01, t02, t03, t04, t05, t06, t07;
  74     t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
  75     t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
  76     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  77     t04 = sqrt(t01 + t02 + t03);
  78     t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
  79     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  80     t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
  81     F += t04;
  82     Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  83     Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
  84     Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
  85     Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  86             2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
  87             2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
  88     Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
  89             2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
  90             2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  91     Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
  92             2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
  93 }
  94
  95 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
  96     double t0 = 0, t1 = 0, t2 = 0;
  97     for (unsigned int i = 0; i < b.size(); i++) {
  98         t0 = static_cast< double >(b[i][0]);
  99         t1 = static_cast< double >(b[i][1]);
 100         t2 = static_cast< double >(b[i][2]);
 101         b[i][0] = static_cast< float >((t2 + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (t1 + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (t0 + x));
 102         b[i][1] = static_cast< float >((t1 + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (t2 + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (t0 + x));
 103         b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
 104     }
 105 }
 106
 107 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
 108     midA[0] = 0;
 109     midA[1] = 0;
 110     midA[2] = 0;
 111
 112     midB[0] = 0;
 113     midB[1] = 0;
 114     midB[2] = 0;
 115
 116     for (unsigned int i = 0; i < pairs.size(); i++) {
 117         midA += a[pairs[i].first];
 118         midB += b[pairs[i].second];
 119         //rvec_inc(midA, a[pairs[i].first]);
 120         //rvec_inc(midB, b[pairs[i].second]);
 121     }
 122     midA[0] /= pairs.size();
 123     midA[1] /= pairs.size();
 124     midA[2] /= pairs.size();
 125
 126     midB[0] /= pairs.size();
 127     midB[1] /= pairs.size();
 128     midB[2] /= pairs.size();
 129 }
 130
 131 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
 132     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 133     RVec ma, mb;
 134     CalcMid(a, b, ma, mb, pairs);
 135     ma -= mb;
 136     //rvec_dec(ma, mb);
 137     double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
 138     for (unsigned int i = 0; i < pairs.size(); i++) {
 139         f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 140                 static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 141     }
 142     if (FtCnst == 0) {
 143         FtCnst = 0.00001;
 144     }
 145     if (f1 == 0) {
 146         return;
 147     } else {
 148         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
 149             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
 150             for (unsigned int i = 0; i < pairs.size(); i++) {
 151                 searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
 152                                static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 153                                static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
 154                                static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 155             }
 156             while (f2 != f1) {
 157                 f2 = 0;
 158                 for (unsigned int i = 0; i < pairs.size(); i++) {
 159                     f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 160                             static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
 161                             static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
 162                 }
 163                 if (f2 < f1) {
 164                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 165                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 166                     break;
 167                 } else {
 168                     l *= 0.50;
 169                     if (l == DBL_MIN) { //DBL_TRUE_MIN
 170                         break;
 171                     }
 172                 }
 173             }
 174         }
 175         ApplyFit(b, x, y, z, A, B, C);
 176     }
 177 }
 178
 179 class domainsType {
 180     public:
 181
 182         ~domainsType() {}
 183
 184         domainsType() {}
 185
 186         domainsType(std::vector< int > index, const std::vector< RVec > reference,
 187                     const int windowSize, const int domainMinimumSize,
 188                     const int domainSearchAlgorythm, const int timeStepBetweenWindowStarts,
 189                     const unsigned int sliceNum, const double epsilon, const double delta,
 190                     const std::string outPutFileName) {
 191             setDefaults(index, reference, windowSize, domainMinimumSize, domainSearchAlgorythm, timeStepBetweenWindowStarts, sliceNum, epsilon, delta, outPutFileName);
 192         }
 193
 194         // set numeric values to "-1" / string value to "" - if you want default settings
 195         void setDefaults(std::vector< int > index, const std::vector< RVec > reference,
 196                         const int windowSize, const int domainMinimumSize,
 197                         const int domainSearchAlgorythm, const int timeStepBetweenWindows,
 198                         const unsigned int sliceNum, const double epsilon, const double delta,
 199                         const std::string outPutFileName) {
 200             graph.resize(1);
 201             graph.front().resize(sliceNum);
 202             for (unsigned i = 0; i < sliceNum; i++) {
 203                 setGraph(graph.front()[i], reference);
 204             }
 205             structIndex = index;
 206             ref = reference;
 207             if (epsilon != -1) { eps = epsilon; }
 208             if (delta != -1) { dlt = delta; }
 209             if (windowSize != -1) { window = windowSize; }
 210             if (domainMinimumSize != -1) { dms = domainMinimumSize; }
 211             if (domainSearchAlgorythm != -1) { dsa = domainSearchAlgorythm; }
 212             if (timeStepBetweenWindows != -1) { ts = timeStepBetweenWindows; }
 213             if (outPutFileName != "") { outPut = outPutFileName; }
 214         }
 215
 216         void update(const std::vector< std::vector< RVec > > frame, const int frameNumber) {
 217             if (frameNumber % ts == 0) {
 218                 graph.resize(graph.size() + 1);
 219                 graph.back().resize(graph.front().size());
 220                 for (unsigned i = 0; i < graph.front().size(); i++) {
 221                     setGraph(graph.back()[i], ref);
 222                 }
 223             }
 224             for (unsigned int i = 0; i < graph.size(); i++) {
 225                 #pragma omp parallel for
 226                 for (unsigned long j = 0; j < graph.front().size(); j++) {
 227                     for (unsigned int k1 = 0; k1 < graph[i][j].size(); k1++) {
 228                         for (unsigned int k2 = k1; k2 < graph[i][j].size(); k2++) {
 229                             if ((frame[j][k1] - frame[j][k2] - graph[i][j][k1][k2].radius).norm() <= static_cast< float >(eps)) {
 230                                 if (k1 == k2) {
 231                                     graph[i][j][k1][k2].num++;
 232                                 } else {
 233                                     graph[i][j][k1][k2].num++;
 234                                     graph[i][j][k2][k1].num++;
 235                                 }
 236                             }
 237                         }
 238                     }
 239                 }
 240                 #pragma omp barrier
 241             }
 242             if (graph.size() > 0) {
 243                 updateCount++;
 244             }
 245             if (updateCount == window) {
 246                 getDomains();
 247                 print(frameNumber);
 248             }
 249         }
 250
 251     private:
 252         struct triple {
 253             int num;
 254             RVec radius;
 255             bool check;
 256         };
 257         std::vector< RVec >                                                     ref;                            // must be presented
 258         std::vector< int >                                                      structIndex;                    // must be presented
 259         std::vector< std::vector< std::vector< std::vector< triple > > > >      graph;
 260         std::vector< std::vector< unsigned int > >                              domains;
 261         std::vector< std::vector< int > >                                       domsizes;
 262         int                                                                     window      = 1000;             // selectable
 263         int                                                                     updateCount = 0;
 264         int                                                                     dms         = 4;                // selectable
 265         int                                                                     dsa         = 0;                // selectable
 266         double                                                                  eps         = 0.2;              // selectable
 267         double                                                                  dlt         = 0.95;             // selectable
 268         int                                                                     ts          = window / 10;      // selectable
 269         std::string                                                             outPut      = "default.ndx";    // selectable
 270
 271         void setGraph(std::vector< std::vector< triple > > &smallGraph, const std::vector< RVec > reference) {
 272             smallGraph.resize(reference.size());
 273             for (unsigned i = 0; i < reference.size(); i++) {
 274                 smallGraph[i].resize(reference.size());
 275                 for (unsigned j = 0; j < reference.size(); j++) {
 276                     smallGraph[i][j].num = 0;
 277                     smallGraph[i][j].radius = reference[i] - reference[j];
 278                     smallGraph[i][j].check = false;
 279                 }
 280             }
 281         }
 282
 283         void deleteDomainFromGraph(std::vector< unsigned int > domain) {
 284             for (unsigned int i = 0; i < graph.front().size(); i++) {
 285                 for (unsigned int j = 0; j < domain.size(); j++) {
 286                     for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 287                         graph.front()[i][domain[j]][k].check = false;
 288                         graph.front()[i][k][domain[j]].check = false;
 289                     }
 290                 }
 291             }
 292         }
 293
 294         bool searchDomainSizes() {
 295             // it goes infinite without these 2 lines, curious even if it was resized before
 296             // resizing inside vectors with 0 is not enough
 297             domsizes.resize(0);
 298             domsizes.resize(graph.front().size());
 299             bool flag = false;
 300             for (unsigned int i = 0; i < graph.front().size(); i++) {
 301                 domsizes[i].resize(graph.front()[i].size(), 0);
 302                 for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
 303                     for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 304                         if (graph.front()[i][j][k].check) {
 305                             domsizes[i][j]++;
 306                         }
 307                         if (!flag) {
 308                             if (domsizes[i][j] >= dms) {
 309                                 flag = true;
 310                             }
 311                         }
 312                     }
 313                 }
 314             }
 315             return flag;
 316         }
 317
 318         void getDomains() {
 319             for (unsigned int i = 0; i < graph.front().size(); i++) {
 320                 for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
 321                     for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 322                         if (graph.front()[i][j][k].num >= window * dlt) {
 323                             graph.front()[i][j][k].check = true;
 324                         }
 325                     }
 326                 }
 327             }
 328             domains.resize(0);
 329             while (searchDomainSizes()) {
 330                 unsigned int t1 = 0, t2 = 0;
 331                 for (unsigned int i = 0; i < domsizes.size(); i++) {
 332                     for (unsigned int j = 0; j < domsizes[i].size(); j++) {
 333                         if ((dsa == 0) && (domsizes[i][j] > domsizes[t1][t2])) {
 334                             // подумать как понять какой домен лучше при одинаковом размере
 335                             t1 = i;
 336                             t2 = j;
 337                         }
 338                         if ((dsa == 1) && ((domsizes[i][j] < domsizes[t1][t2]) || ((domsizes[i][j] >= dms) && (domsizes[t1][t2] < dms)))) {
 339                             t1 = i;
 340                             t2 = j;
 341                         }
 342                     }
 343                 }
 344                 if (domsizes[t1][t2] >= dms) {
 345                     domains.resize(domains.size() + 1);
 346                     for (unsigned int i = 0; i < graph.front()[t1][t2].size(); i++) {
 347                         if (graph.front()[t1][t2][i].check) {
 348                             domains.back().push_back(i);
 349                         }
 350                     }
 351                     deleteDomainFromGraph(domains.back());
 352                 } else {
 353                     break;
 354                 }
 355             }
 356             graph.erase(graph.begin());
 357             updateCount = std::max(0, window - ts);
 358         }
 359
 360         void print(int currentFrame) {
 361             if (domains.size() == 0) {
 362                 return;
 363             }
 364             FILE                *ndxFile, *slFile;
 365             ndxFile = std::fopen(outPut.c_str(), "a+");
 366             slFile = std::fopen(("SelectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+");
 367             int writeCount;
 368             for (unsigned int i = 0; i < domains.size(); i++) {
 369                     std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", currentFrame - window, i + 1, dms, eps, dlt);
 370                     std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', currentFrame - window, i + 1, dms, eps, dlt, '"');
 371                     writeCount = 0;
 372                     for (unsigned int j = 0; j < domains[i].size(); j++) {
 373                         writeCount++;
 374                         if (writeCount > 20) {
 375                             writeCount -= 20;
 376                             std::fprintf(ndxFile, "\n");
 377                         }
 378                         std::fprintf(ndxFile, "%5d ", structIndex[domains[i][j]] + 1);
 379                         std::printf("%5d ", structIndex[domains[i][j]] + 1);
 380                     }
 381                     std::fprintf(ndxFile,"\n\n");
 382                     std::printf("\n\n");
 383             }
 384             std::fprintf(ndxFile,"\n");
 385             std::fclose(ndxFile);
 386             std::fclose(slFile);
 387         }
 388 };
 389
 390 struct node {
 391     short int n;
 392     RVec r;
 393     bool yep;
 394 };
 395
 396 void make_graph(unsigned long mgwi_natoms, std::vector < RVec > mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
 397 {
 398     mgwi_graph.resize(mgwi_natoms);
 399     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 400         mgwi_graph[i].resize(mgwi_natoms);
 401     }
 402     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 403         for (unsigned int j = 0; j < mgwi_natoms; j++) {
 404             rvec_sub(mgwi_x[i].as_vec(), mgwi_x[j].as_vec(), mgwi_graph[i][j].r.as_vec());
 405             mgwi_graph[i][j].n = 0;
 406         }
 407     }
 408 }
 409
 410 void update_graph(std::vector< std::vector< node > > &ugwi_graph, std::vector < RVec > ugwi_x, /*long*/ double ugwi_epsi) {
 411     RVec ugwi_temp;
 412     for (unsigned int i = 0; i < ugwi_graph.size(); i++) {
 413         for (unsigned int j = i; j < ugwi_graph.size(); j++) {
 414             rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), ugwi_temp.as_vec());
 415             rvec_dec(ugwi_temp.as_vec(), ugwi_graph[i][j].r.as_vec());
 416             if (static_cast< double >(norm(ugwi_temp)) <= ugwi_epsi) {
 417                 if (i == j) {
 418                     ugwi_graph[i][j].n++;
 419                 }
 420                 else {
 421                     ugwi_graph[i][j].n++;
 422                     ugwi_graph[j][i].n++;
 423                 }
 424             }
 425         }
 426     }
 427 }
 428
 429 void check_domains(double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
 430     for (unsigned int k = 0; k < cd_graph.size(); k++) {
 431         for (unsigned int i = 0; i < cd_graph[1].size(); i++) {
 432             for (unsigned int j = 0; j < cd_graph[1].size(); j++) {
 433                 if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
 434                     cd_graph[k][i][j].yep = true;
 435                 }
 436                 else {
 437                     cd_graph[k][i][j].yep = false;
 438                 }
 439             }
 440         }
 441     }
 442 }
 443
 444 void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
 445     fds_domsizes.resize(fds_graph.size());
 446     for (unsigned int i = 0; i < fds_graph.size(); i++) {
 447         fds_domsizes[i].resize(fds_graph[1].size(), 0);
 448         for (unsigned int j = 0; j < fds_graph[1].size(); j++) {
 449             for (unsigned int k = 0; k < fds_graph[1].size(); k++) {
 450                 if (fds_graph[i][j][k].yep) {
 451                     fds_domsizes[i][j]++;
 452                 }
 453             }
 454         }
 455     }
 456 }
 457
 458 void get_maxsized_domain(std::vector< unsigned int > &gmd_max_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes) {
 459     unsigned int gmd_number1 = 0, gmd_number2 = 0;
 460     for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
 461         for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
 462             if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
 463                 gmd_number1 = i;
 464                 gmd_number2 = j;
 465             }
 466         }
 467     }
 468     gmd_max_d.resize(0);
 469     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 470         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 471             gmd_max_d.push_back(i);
 472         }
 473     }
 474 }
 475
 476 void get_min_domain(std::vector< unsigned int > &gmd_min_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes, int gmd_min_dom_size) {
 477     unsigned int gmd_number1 = 0, gmd_number2 = 0;
 478     for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
 479         for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
 480             if (gmd_domsizes[gmd_number1][gmd_number2] < gmd_min_dom_size && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 481                 gmd_number1 = i;
 482                 gmd_number2 = j;
 483             }
 484             if (gmd_domsizes[i][j] <= gmd_domsizes[gmd_number1][gmd_number2] && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 485                 gmd_number1 = i;
 486                 gmd_number2 = j;
 487             }
 488         }
 489     }
 490     gmd_min_d.resize(0);
 491     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 492         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 493             gmd_min_d.push_back(i);
 494         }
 495     }
 496 }
 497
 498 void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< unsigned int > ddf_domain) {
 499     for (unsigned int i = 0; i < ddf_domain.size(); i++) {
 500         for (unsigned int k = 0; k < ddf_graph.size(); k++) {
 501             for (unsigned int j = 0; j < ddf_graph[1].size(); j++) {
 502                 if (ddf_graph[k][ddf_domain[i]][j].yep) {
 503                     ddf_graph[k][ddf_domain[i]][j].yep = false;
 504                 }
 505                 if (ddf_graph[k][j][ddf_domain[i]].yep) {
 506                     ddf_graph[k][j][ddf_domain[i]].yep = false;
 507                 }
 508             }
 509         }
 510     }
 511 }
 512
 513 bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
 514     for (unsigned int i = 0; i < cd_domsizes.size(); i++) {
 515         for (unsigned int j = 0; j < cd_domsizes[0].size(); j++) {
 516             if (cd_domsizes[i][j] >= cd_domain_min_size) {
 517                 return true;
 518             }
 519         }
 520     }
 521     return false;
 522 }
 523
 524 void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta, int window, int st_pos) {
 525     FILE                *fpNdx_, *fpNdx2_;
 526     fpNdx_ = std::fopen(fnNdx_.c_str(), "a+");
 527     fpNdx2_ = std::fopen(("SelectionList-" + fnNdx_.substr(0, fnNdx_.size() - 4)).c_str(), "a+");
 528     int write_count;
 529     for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) {
 530             std::fprintf(fpNdx_, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta);
 531             std::fprintf(fpNdx2_, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta, '"');
 532             write_count = 0;
 533             for (unsigned int j = 0; j < pd_domains[i1].size(); j++) {
 534                 write_count++;
 535                 if (write_count > 20) {
 536                     write_count -= 20;
 537                     std::fprintf(fpNdx_, "\n");
 538                 }
 539                 std::fprintf(fpNdx_, "%5d ", index[pd_domains[i1][j]] + 1);
 540                 std::printf("%5d ", index[pd_domains[i1][j]] + 1);
 541             }
 542             std::printf("\n\n");
 543             std::fprintf(fpNdx_,"\n\n");
 544     }
 545     std::fprintf(fpNdx_,"\n");
 546     std::fclose(fpNdx_);
 547     //std::fprintf(fpNdx2_,"\n");
 548     std::fclose(fpNdx2_);
 549 }
 550
 551 /*! \brief
 552  * \ingroup module_trajectoryanalysis
 553  */
 554 class Domains : public TrajectoryAnalysisModule
 555 {
 556     public:
 557
 558         Domains();
 559         virtual ~Domains();
 560
 561         //! Set the options and setting
 562         virtual void initOptions(IOptionsContainer          *options,
 563                                  TrajectoryAnalysisSettings *settings);
 564
 565         //! First routine called by the analysis framework
 566         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 567         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 568                                   const TopologyInformation        &top);
 569
 570         //! Call for each frame of the trajectory
 571         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 572         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 573                                   TrajectoryAnalysisModuleData *pdata);
 574
 575         //! Last routine called by the analysis framework
 576         // virtual void finishAnalysis(t_pbc *pbc);
 577         virtual void finishAnalysis(int nframes);
 578
 579         //! Routine to write output, that is additional over the built-in
 580         virtual void writeOutput();
 581
 582     private:
 583
 584         std::string                                                             fnNdx_;
 585
 586         std::vector< std::vector< std::vector< std::vector< node > > > >        graph;
 587
 588         domainsType                                                             testSubject;
 589
 590
 591         std::vector< std::vector< unsigned int > >                              domains;
 592         std::vector< std::vector< int > >                                       domsizes;
 593
 594         std::vector< int >                                                      index;
 595         std::vector< RVec >                                                     trajectory;
 596         std::vector< RVec >                                                     reference;
 597         Selection                                                               selec;
 598         int                                                                     frames              = 0;
 599         int                                                                     window              = -1; // selectable
 600         int                                                                     domain_min_size     = -1; // selectable
 601
 602         std::vector< std::vector< std::pair< unsigned int, unsigned int > > >   w_rls2;
 603         unsigned int                                                            bone;
 604         double                                                                  delta               = -1; // selectable
 605         double                                                                  epsi                = -1; // selectable
 606         int                                                                     twStep              = -1;
 607
 608         int                                                                     DomainSearchingAlgorythm = -1; // selectable
 609         // Copy and assign disallowed by base.
 610 };
 611
 612 Domains::Domains(): TrajectoryAnalysisModule()
 613 {
 614 }
 615
 616 Domains::~Domains()
 617 {
 618 }
 619
 620 void
 621 Domains::initOptions(IOptionsContainer          *options,
 622                      TrajectoryAnalysisSettings *settings)
 623 {
 624     static const char *const desc[] = {
 625         "[THISMODULE] to be done"
 626     };
 627     // Add the descriptive text (program help text) to the options
 628     settings->setHelpText(desc);
 629     // Add option for selecting a subset of atoms
 630     options->addOption(SelectionOption("select")
 631                            .store(&selec).required()
 632                            .description("Atoms that are considered as part of the excluded volume"));
 633     // Add option for output file name
 634     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 635                             .store(&fnNdx_).defaultBasename("domains")
 636                             .description("Index file from the domains"));
 637     // Add option for domain min size constant
 638     options->addOption(gmx::IntegerOption("dms")
 639                             .store(&domain_min_size)
 640                             .description("minimum domain size, should be >= 4"));
 641     // Add option for Domain's Searching Algorythm
 642     options->addOption(gmx::IntegerOption("DSA")
 643                             .store(&DomainSearchingAlgorythm)
 644                             .description("Domain's Searching Algorythm: 0 == default (from bigger to smaller) | 1 == (from smaller to bigger)"));
 645     // Add option for epsi constant
 646     options->addOption(DoubleOption("eps")
 647                             .store(&epsi)
 648                             .description("thermal vibrations' constant"));
 649     // Add option for delta constant
 650     options->addOption(DoubleOption("dlt")
 651                             .store(&delta)
 652                             .description("domain membership probability"));
 653     // Add option for domain min size constant
 654     options->addOption(gmx::IntegerOption("wSize")
 655                             .store(&window)
 656                             .description("flowing window to evaluate domains from"));
 657     // Add option for domain min size constant
 658     options->addOption(gmx::IntegerOption("twStep")
 659                             .store(&twStep)
 660                             .description("time step between windows' starting positions"));
 661     // Control input settings
 662     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 663     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 664     settings->setPBC(true);
 665 }
 666
 667 void
 668 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 669                       const TopologyInformation        &top)
 670 {
 671     index.resize(0);
 672     for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
 673         index.push_back(*ai);
 674     }
 675
 676     reference.resize(0);
 677     if (top.hasFullTopology()) {
 678         for (unsigned int i = 0; i < index.size(); i++) {
 679             reference.push_back(top.x().at(index[i]));
 680         }
 681     }
 682     bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
 683
 684     graph.resize(1);
 685     graph.front().resize(bone);
 686
 687     w_rls2.resize(bone);
 688     for (unsigned int i = 0; i < bone; i++) {
 689         w_rls2[i].resize(0);
 690         for (unsigned int j = 0; j < static_cast< unsigned int >(domain_min_size); j++) {
 691             w_rls2[i].push_back(std::make_pair(j + i, j + i));
 692         }
 693         //make_graph(index.size(), reference, graph.front()[i]);
 694     }
 695
 696     trajectory.resize(index.size());
 697
 698     testSubject.setDefaults(index, reference, window, domain_min_size, DomainSearchingAlgorythm, twStep, bone, epsi, delta, fnNdx_);
 699 }
 700
 701 void
 702 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 703 {
 704     /*if (frnr % twStep == 0) {
 705         graph.resize(graph.size() + 1);
 706         graph.back().resize(bone);
 707         for (unsigned int i = 0; i < bone; i++) {
 708             make_graph(index.size(), reference, graph.back()[i]);
 709         }
 710     }
 711
 712     int temp = graph.size() - window / twStep;
 713
 714     // most probably == 0
 715     if (temp > 0) {
 716         domains.resize(0);
 717         std::vector< unsigned int > a;
 718         a.resize(0);
 719         domsizes.resize(0);
 720         check_domains(delta, window, graph[0]);
 721         std::cout << "finding domains' sizes\n";
 722         find_domain_sizes(graph[0], domsizes);
 723         while (check_domsizes(domsizes, domain_min_size)) {
 724             domains.push_back(a);
 725             if (DomainSearchingAlgorythm == 0) {
 726                 get_maxsized_domain(domains.back(), graph[0], domsizes);
 727             } else if (DomainSearchingAlgorythm == 1) {
 728                 get_min_domain(domains.back(), graph[0], domsizes, domain_min_size);
 729             }
 730             std::cout << "new domain: " << domains.back().size() << " atoms\n";
 731             delete_domain_from_graph(graph[0], domains.back());
 732             domsizes.resize(0);
 733             find_domain_sizes(graph[0], domsizes);
 734         }
 735         graph.erase(graph.begin());
 736         std::cout << (frnr + 1) / twStep - (window / twStep)<< " window's domains have been evaluated\n";
 737         print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / twStep - (window / twStep)); // see function for details | numbers from index
 738     }*/
 739
 740     for (unsigned int i = 0; i < index.size(); i++) {
 741         trajectory[i] = fr.x[index[i]];
 742     }
 743     frames++;
 744
 745     /*#pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
 746     for (unsigned int j = 0; j < bone; j++) {
 747         std::vector < RVec > temp = trajectory;
 748         MyFitNew(reference, temp, w_rls2[j], 0);
 749         for (unsigned int i = 0; i < graph.size(); i++) {
 750             update_graph(graph[i][j], temp, epsi);
 751         }
 752     }
 753     #pragma omp barrier*/
 754
 755     std::vector< std::vector< RVec > > tsTemp;
 756     tsTemp.resize(bone, trajectory);
 757
 758     #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2) firstprivate(trajectory, reference)
 759     for (unsigned int i = 0; i < bone; i++) {
 760         MyFitNew(reference, tsTemp[i], w_rls2[i], 0);
 761     }
 762     #pragma omp barrier
 763     testSubject.update(tsTemp, frnr);
 764     //std::cout << "frame: " << frnr << " analyzed\n";
 765 }
 766
 767 void
 768 Domains::finishAnalysis(int nframes)
 769 {
 770
 771 }
 772
 773 void
 774 Domains::writeOutput()
 775 {
 776     std::cout << "\n END \n";
 777 }
 778
 779 /*! \brief
 780  * The main function for the analysis template.
 781  */
 782 int
 783 main(int argc, char *argv[])
 784 {
 785     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 786 }