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 const int tau_jump = 10;
  67
  68 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) {
  69     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) +
  70                     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) +
  71                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  72 }
  73
  74 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) {
  75     double t01, t02, t03, t04, t05, t06, t07;
  76     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);
  77     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);
  78     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  79     t04 = sqrt(t01 + t02 + t03);
  80     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);
  81     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  82     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);
  83     F += t04;
  84     Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  85     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);
  86     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);
  87     Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  88             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 +
  89             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);
  90     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 +
  91             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 -
  92             2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  93     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 -
  94             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);
  95 }
  96
  97 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
  98     double t0 = 0, t1 = 0, t2 = 0;
  99     for (unsigned int i = 0; i < b.size(); i++) {
 100         t0 = static_cast< double >(b[i][0]);
 101         t1 = static_cast< double >(b[i][1]);
 102         t2 = static_cast< double >(b[i][2]);
 103         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));
 104         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));
 105         b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
 106     }
 107 }
 108
 109 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
 110     midA[0] = 0;
 111     midA[1] = 0;
 112     midA[2] = 0;
 113
 114     midB[0] = 0;
 115     midB[1] = 0;
 116     midB[2] = 0;
 117
 118     for (unsigned int i = 0; i < pairs.size(); i++) {
 119         midA += a[pairs[i].first];
 120         midB += b[pairs[i].second];
 121         //rvec_inc(midA, a[pairs[i].first]);
 122         //rvec_inc(midB, b[pairs[i].second]);
 123     }
 124     midA[0] /= pairs.size();
 125     midA[1] /= pairs.size();
 126     midA[2] /= pairs.size();
 127
 128     midB[0] /= pairs.size();
 129     midB[1] /= pairs.size();
 130     midB[2] /= pairs.size();
 131 }
 132
 133 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
 134     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 135     RVec ma, mb;
 136     CalcMid(a, b, ma, mb, pairs);
 137     ma -= mb;
 138     //rvec_dec(ma, mb);
 139     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;
 140     for (unsigned int i = 0; i < pairs.size(); i++) {
 141         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]),
 142                 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);
 143     }
 144     if (FtCnst == 0) {
 145         FtCnst = 0.00001;
 146     }
 147     if (f1 == 0) {
 148         return;
 149     } else {
 150         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
 151             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
 152             for (unsigned int i = 0; i < pairs.size(); i++) {
 153                 searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
 154                                static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 155                                static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
 156                                static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 157             }
 158             while (f2 != f1) {
 159                 f2 = 0;
 160                 for (unsigned int i = 0; i < pairs.size(); i++) {
 161                     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]),
 162                             static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
 163                             static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
 164                 }
 165                 if (f2 < f1) {
 166                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 167                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 168                     break;
 169                 } else {
 170                     l *= 0.50;
 171                     if (l == DBL_MIN) { //DBL_TRUE_MIN
 172                         break;
 173                     }
 174                 }
 175             }
 176         }
 177         ApplyFit(b, x, y, z, A, B, C);
 178     }
 179 }
 180
 181 class domainsType {
 182     public:
 183
 184         ~domainsType() {}
 185
 186         domainsType() {}
 187
 188         domainsType(const unsigned int sliceNum, const std::vector< RVec > reference) {
 189             setDefault(sliceNum, reference);
 190         }
 191
 192         void update(const std::vector< std::vector< RVec > > frame, const float epsilon, const int frameNumber) {
 193             if (updateCount == window)
 194             if ((frameNumber + 1) % ts == 0) {
 195                 graph.resize(graph.size() + 1);
 196                 graph.back().resize(graph.front().size());
 197                 for (unsigned i = 0; i < graph.front().size(); i++) {
 198                     setGraph(graph.back()[i], ref);
 199                 }
 200             }
 201             for (unsigned int i = 0; i < graph.size(); i++) {
 202                 for (unsigned int j = 0; j < graph.front().size(); j++) {
 203                     for (unsigned int k1 = 0; k1 < graph[i][j].size(); k1++) {
 204                         for (unsigned int k2 = k1; k2 < graph[i][j].size(); k2++) {
 205                             if ((frame[j][k1] - frame[j][k2] - graph[i][j][k1][k2].radius).norm() <= epsilon) {
 206                                 if (k1 == k2) {
 207                                     graph[i][j][k1][k2].num++;
 208                                 } else {
 209                                     graph[i][j][k1][k2].num++;
 210                                     graph[i][j][k2][k1].num++;
 211                                 }
 212                             }
 213                         }
 214                     }
 215                 }
 216             }
 217             updateCount++;
 218         }
 219
 220         void getDomains(const float delta) {
 221             if (updateCount != window) {
 222                 return;
 223             } else {
 224                 //присмотреться к этому моменту
 225                 updateCount -= ts;
 226                 for (unsigned int i = 0; i < graph.front().size(); i++) {
 227                     for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
 228                         for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 229                             if (graph.front()[i][j][k].num >= window * delta) {
 230                                 graph.front()[i][j][k].check = true;
 231                             }
 232                         }
 233                     }
 234                 }
 235                 domains.resize(0);
 236                 while (checkDomainSizes()) {
 237                     domsizes.resize(graph.front().size());
 238                     for (unsigned int i = 0; i < graph.front().size(); i++) {
 239                         domsizes[i].resize(graph.front()[i].size(), 0);
 240                         for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
 241                             for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 242                                 if (graph.front()[i][j][k].check) {
 243                                     domsizes[i][j]++;
 244                                 }
 245                             }
 246                         }
 247                     }
 248                     unsigned int t1 = 0, t2 = 0;
 249                     for (unsigned int i = 0; i < domsizes.size(); i++) {
 250                         for (unsigned int j = 0; j < domsizes[i].size(); j++) {
 251                             if ((dsa == 0) && (domsizes[i][j] > domsizes[t1][t2])) {
 252                                 // подумать как понять какой домен лучше при одинаковом размере
 253                                 t1 = i;
 254                                 t2 = j;
 255                             }
 256                             if ((dsa == 1) && ((domsizes[i][j] < domsizes[t1][t2]) || ((domsizes[i][j] >= dms) && (domsizes[t1][t2] < dms)))) {
 257                                 t1 = i;
 258                                 t2 = j;
 259                             }
 260                         }
 261                     }
 262                     domains.resize(domains.size() + 1);
 263                     for (unsigned int i = 0; i < graph.front()[t1][t2].size(); i++) {
 264                         if (graph.front()[t1][t2][i].check) {
 265                             domains.back().push_back(i);
 266                         }
 267                     }
 268                     deleteDomainFromGraph(domains.back());
 269                 }
 270                 if (graph.size() > static_cast< unsigned int >(window / ts)) {
 271                     graph.erase(graph.begin());
 272                 }
 273             }
 274         }
 275
 276         void setParams(const int windowSize, const int domainMinimumSize, const int domainSearchAlgorythm, const int timeStepBetweenWindows) {
 277             window = windowSize;
 278             dms = domainMinimumSize;
 279             dsa = domainSearchAlgorythm;
 280             ts = timeStepBetweenWindows;
 281         }
 282
 283         void setDefault(const unsigned int sliceNum, const std::vector< RVec > reference) {
 284             graph.resize(1);
 285             graph.front().resize(sliceNum);
 286             for (unsigned i = 0; i < sliceNum; i++) {
 287                 setGraph(graph.front()[i], reference);
 288             }
 289             ref = reference;
 290         }
 291
 292         void print(std::vector< int > index, std::string fileName, double epsilon, double delta, int startingPos) {
 293             if (domains.size() == 0) {
 294                 return;
 295             }
 296             FILE                *ndxFile, *slFile;
 297             ndxFile = std::fopen(fileName.c_str(), "a+");
 298             slFile = std::fopen(("SelectionList-" + fileName.substr(0, fileName.size() - 4)).c_str(), "a+");
 299             int writeCount;
 300             for (unsigned int i = 0; i < domains.size(); i++) {
 301                     std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(startingPos) * ts, i + 1, dms, epsilon, delta);
 302                     std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(startingPos) * ts, i + 1, dms, epsilon, delta, '"');
 303                     writeCount = 0;
 304                     for (unsigned int j = 0; j < domains[i].size(); j++) {
 305                         writeCount++;
 306                         if (writeCount > 20) {
 307                             writeCount -= 20;
 308                             std::fprintf(ndxFile, "\n");
 309                         }
 310                         std::fprintf(ndxFile, "%5d ", index[domains[i][j]] + 1);
 311                         std::printf("%5d ", index[domains[i][j]] + 1);
 312                     }
 313                     std::fprintf(ndxFile,"\n\n");
 314                     std::printf("\n\n");
 315             }
 316             std::fprintf(ndxFile,"\n");
 317             std::fclose(ndxFile);
 318             std::fclose(slFile);
 319         }
 320
 321     private:
 322         struct triple {
 323             int num;
 324             RVec radius;
 325             bool check;
 326         };
 327         std::vector< RVec >                                                     ref;
 328         std::vector< std::vector< std::vector< std::vector< triple > > > >      graph;
 329         std::vector< std::vector< unsigned int > >                              domains;
 330         std::vector< std::vector< int > >                                       domsizes;
 331         int                                                                     window      = 1000; // selectable
 332         int                                                                     updateCount = 0;
 333         int                                                                     dms         = 4;    // selectable
 334         int                                                                     dsa         = 0;    // selectable
 335         int                                                                     ts          = window / 10;
 336
 337         void setGraph(std::vector< std::vector< triple > > &smallGraph, const std::vector< RVec > reference) {
 338             smallGraph.resize(reference.size());
 339             for (unsigned i = 0; i < reference.size(); i++) {
 340                 smallGraph[i].resize(reference.size());
 341                 for (unsigned j = 0; j < reference.size(); j++) {
 342                     smallGraph[i][j].num = 0;
 343                     smallGraph[i][j].radius = reference[i] - reference[j];
 344                     smallGraph[i][j].check = false;
 345                 }
 346             }
 347         }
 348
 349         void deleteDomainFromGraph(std::vector< unsigned int > domain) {
 350             for (unsigned int i = 0; i < graph.front().size(); i++) {
 351                 for (unsigned int j = 0; j < domain.size(); j++) {
 352                     for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
 353                         graph.front()[i][domain[j]][k].check = false;
 354                         graph.front()[i][k][domain[j]].check = false;
 355                     }
 356                 }
 357             }
 358         }
 359
 360         bool checkDomainSizes() {
 361             for (unsigned int i = 0; i < domsizes.size(); i++) {
 362                 for (unsigned int j = 0; j < domsizes[i].size(); j++) {
 363                     if (domsizes[i][j] >= dms) {
 364                         return true;
 365                     }
 366                 }
 367             }
 368             return false;
 369         }
 370 };
 371
 372 struct node {
 373     short int n;
 374     RVec r;
 375     bool yep;
 376 };
 377
 378 void make_graph(unsigned long mgwi_natoms, std::vector < RVec > mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
 379 {
 380     mgwi_graph.resize(mgwi_natoms);
 381     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 382         mgwi_graph[i].resize(mgwi_natoms);
 383     }
 384     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 385         for (unsigned int j = 0; j < mgwi_natoms; j++) {
 386             rvec_sub(mgwi_x[i].as_vec(), mgwi_x[j].as_vec(), mgwi_graph[i][j].r.as_vec());
 387             mgwi_graph[i][j].n = 0;
 388         }
 389     }
 390 }
 391
 392 void update_graph(std::vector< std::vector< node > > &ugwi_graph, std::vector < RVec > ugwi_x, /*long*/ double ugwi_epsi) {
 393     RVec ugwi_temp;
 394     for (unsigned int i = 0; i < ugwi_graph.size(); i++) {
 395         for (unsigned int j = i; j < ugwi_graph.size(); j++) {
 396             rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), ugwi_temp.as_vec());
 397             rvec_dec(ugwi_temp.as_vec(), ugwi_graph[i][j].r.as_vec());
 398             if (static_cast< double >(norm(ugwi_temp)) <= ugwi_epsi) {
 399                 if (i == j) {
 400                     ugwi_graph[i][j].n++;
 401                 }
 402                 else {
 403                     ugwi_graph[i][j].n++;
 404                     ugwi_graph[j][i].n++;
 405                 }
 406             }
 407         }
 408     }
 409 }
 410
 411 void check_domains(double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
 412     for (unsigned int k = 0; k < cd_graph.size(); k++) {
 413         for (unsigned int i = 0; i < cd_graph[1].size(); i++) {
 414             for (unsigned int j = 0; j < cd_graph[1].size(); j++) {
 415                 if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
 416                     cd_graph[k][i][j].yep = true;
 417                 }
 418                 else {
 419                     cd_graph[k][i][j].yep = false;
 420                 }
 421             }
 422         }
 423     }
 424 }
 425
 426 void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
 427     fds_domsizes.resize(fds_graph.size());
 428     for (unsigned int i = 0; i < fds_graph.size(); i++) {
 429         fds_domsizes[i].resize(fds_graph[1].size(), 0);
 430         for (unsigned int j = 0; j < fds_graph[1].size(); j++) {
 431             for (unsigned int k = 0; k < fds_graph[1].size(); k++) {
 432                 if (fds_graph[i][j][k].yep) {
 433                     fds_domsizes[i][j]++;
 434                 }
 435             }
 436         }
 437     }
 438 }
 439
 440 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) {
 441     unsigned int gmd_number1 = 0, gmd_number2 = 0;
 442     for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
 443         for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
 444             if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
 445                 gmd_number1 = i;
 446                 gmd_number2 = j;
 447             }
 448         }
 449     }
 450     gmd_max_d.resize(0);
 451     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 452         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 453             gmd_max_d.push_back(i);
 454         }
 455     }
 456 }
 457
 458 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) {
 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[gmd_number1][gmd_number2] < gmd_min_dom_size && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 463                 gmd_number1 = i;
 464                 gmd_number2 = j;
 465             }
 466             if (gmd_domsizes[i][j] <= gmd_domsizes[gmd_number1][gmd_number2] && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 467                 gmd_number1 = i;
 468                 gmd_number2 = j;
 469             }
 470         }
 471     }
 472     gmd_min_d.resize(0);
 473     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 474         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 475             gmd_min_d.push_back(i);
 476         }
 477     }
 478 }
 479
 480 void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< unsigned int > ddf_domain) {
 481     for (unsigned int i = 0; i < ddf_domain.size(); i++) {
 482         for (unsigned int k = 0; k < ddf_graph.size(); k++) {
 483             for (unsigned int j = 0; j < ddf_graph[1].size(); j++) {
 484                 if (ddf_graph[k][ddf_domain[i]][j].yep) {
 485                     ddf_graph[k][ddf_domain[i]][j].yep = false;
 486                 }
 487                 if (ddf_graph[k][j][ddf_domain[i]].yep) {
 488                     ddf_graph[k][j][ddf_domain[i]].yep = false;
 489                 }
 490             }
 491         }
 492     }
 493 }
 494
 495 bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
 496     for (unsigned int i = 0; i < cd_domsizes.size(); i++) {
 497         for (unsigned int j = 0; j < cd_domsizes[0].size(); j++) {
 498             if (cd_domsizes[i][j] >= cd_domain_min_size) {
 499                 return true;
 500             }
 501         }
 502     }
 503     return false;
 504 }
 505
 506 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) {
 507     FILE                *fpNdx_, *fpNdx2_;
 508     fpNdx_ = std::fopen(fnNdx_.c_str(), "a+");
 509     fpNdx2_ = std::fopen(("SelectionList-" + fnNdx_.substr(0, fnNdx_.size() - 4)).c_str(), "a+");
 510     int write_count;
 511     for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) {
 512             std::fprintf(fpNdx_, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(st_pos) * window / tau_jump, i1 + 1, dms, epsi, delta);
 513             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 / tau_jump, i1 + 1, dms, epsi, delta, '"');
 514             write_count = 0;
 515             for (unsigned int j = 0; j < pd_domains[i1].size(); j++) {
 516                 write_count++;
 517                 if (write_count > 20) {
 518                     write_count -= 20;
 519                     std::fprintf(fpNdx_, "\n");
 520                 }
 521                 std::fprintf(fpNdx_, "%5d ", index[pd_domains[i1][j]] + 1);
 522                 std::printf("%5d ", index[pd_domains[i1][j]] + 1);
 523             }
 524             std::printf("\n\n");
 525             std::fprintf(fpNdx_,"\n\n");
 526     }
 527     std::fprintf(fpNdx_,"\n");
 528     std::fclose(fpNdx_);
 529     //std::fprintf(fpNdx2_,"\n");
 530     std::fclose(fpNdx2_);
 531 }
 532
 533 /*! \brief
 534  * \ingroup module_trajectoryanalysis
 535  */
 536 class Domains : public TrajectoryAnalysisModule
 537 {
 538     public:
 539
 540         Domains();
 541         virtual ~Domains();
 542
 543         //! Set the options and setting
 544         virtual void initOptions(IOptionsContainer          *options,
 545                                  TrajectoryAnalysisSettings *settings);
 546
 547         //! First routine called by the analysis framework
 548         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 549         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 550                                   const TopologyInformation        &top);
 551
 552         //! Call for each frame of the trajectory
 553         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 554         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 555                                   TrajectoryAnalysisModuleData *pdata);
 556
 557         //! Last routine called by the analysis framework
 558         // virtual void finishAnalysis(t_pbc *pbc);
 559         virtual void finishAnalysis(int nframes);
 560
 561         //! Routine to write output, that is additional over the built-in
 562         virtual void writeOutput();
 563
 564     private:
 565
 566         std::string                                                             fnNdx_;
 567
 568         std::vector< std::vector< std::vector< std::vector< node > > > >        graph;
 569
 570         domainsType                                                             testSubject;
 571
 572
 573         std::vector< std::vector< unsigned int > >                              domains;
 574         std::vector< std::vector< int > >                                       domsizes;
 575
 576         std::vector< int >                                                      index;
 577         std::vector< RVec >                                                     trajectory;
 578         std::vector< RVec >                                                     reference;
 579         Selection                                                               selec;
 580         int                                                                     frames              = 0;
 581         int                                                                     window              = 1000; // selectable
 582         int                                                                     domain_min_size     = 4; // selectable
 583
 584         std::vector< std::vector< std::pair< unsigned int, unsigned int > > >   w_rls2;
 585         unsigned int                                                            bone;
 586         double                                                                  delta               = 0.95; // selectable
 587         double                                                                  epsi                = 0.10; // selectable
 588
 589         int                                                                     DomainSearchingAlgorythm = 0; // selectable
 590         // Copy and assign disallowed by base.
 591 };
 592
 593 Domains::Domains(): TrajectoryAnalysisModule()
 594 {
 595 }
 596
 597 Domains::~Domains()
 598 {
 599 }
 600
 601 void
 602 Domains::initOptions(IOptionsContainer          *options,
 603                      TrajectoryAnalysisSettings *settings)
 604 {
 605     static const char *const desc[] = {
 606         "[THISMODULE] to be done"
 607     };
 608     // Add the descriptive text (program help text) to the options
 609     settings->setHelpText(desc);
 610     // Add option for selecting a subset of atoms
 611     options->addOption(SelectionOption("select")
 612                            .store(&selec).required()
 613                            .description("Atoms that are considered as part of the excluded volume"));
 614     // Add option for output file name
 615     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 616                             .store(&fnNdx_).defaultBasename("domains")
 617                             .description("Index file from the domains"));
 618     // Add option for domain min size constant
 619     options->addOption(gmx::IntegerOption("dms")
 620                             .store(&domain_min_size)
 621                             .description("minimum domain size, should be >= 4"));
 622     // Add option for Domain's Searching Algorythm
 623     options->addOption(gmx::IntegerOption("DSA")
 624                             .store(&DomainSearchingAlgorythm)
 625                             .description("Domain's Searching Algorythm: 0 == default (from bigger to smaller) | 1 == (from smaller to bigger)"));
 626     // Add option for epsi constant
 627     options->addOption(DoubleOption("epsilon")
 628                             .store(&epsi)
 629                             .description("thermal vibrations' constant"));
 630     // Add option for delta constant
 631     options->addOption(DoubleOption("delta")
 632                             .store(&delta)
 633                             .description("domain membership probability"));
 634     // Control input settings
 635     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 636     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 637     settings->setPBC(true);
 638 }
 639
 640 void
 641 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 642                       const TopologyInformation        &top)
 643 {
 644     index.resize(0);
 645     for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
 646         index.push_back(*ai);
 647     }
 648
 649     reference.resize(0);
 650     if (top.hasFullTopology()) {
 651         for (unsigned int i = 0; i < index.size(); i++) {
 652             reference.push_back(top.x().at(index[i]));
 653         }
 654     }
 655     bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
 656
 657     graph.resize(1);
 658     graph.front().resize(bone);
 659
 660     w_rls2.resize(bone);
 661     for (unsigned int i = 0; i < bone; i++) {
 662         w_rls2[i].resize(0);
 663         for (unsigned int j = 0; j < static_cast< unsigned int >(domain_min_size); j++) {
 664             w_rls2[i].push_back(std::make_pair(j + i, j + i));
 665         }
 666         make_graph(index.size(), reference, graph.front()[i]);
 667     }
 668
 669     trajectory.resize(index.size());
 670
 671     //testSubject.setDefault(bone, reference);
 672     //testSubject.setParams(window, domain_min_size, DomainSearchingAlgorythm, window / tau_jump);
 673 }
 674
 675 void
 676 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 677 {
 678     if ((frnr + 1) % (window / tau_jump) == 0) {
 679         graph.resize(graph.size() + 1);
 680         graph.back().resize(bone);
 681         for (unsigned int i = 0; i < bone; i++) {
 682             make_graph(index.size(), reference, graph.back()[i]);
 683         }
 684     }
 685
 686     int temp = graph.size() - tau_jump;
 687
 688     if (temp > 0) {
 689         domains.resize(0);
 690         std::vector< unsigned int > a;
 691         a.resize(0);
 692         domsizes.resize(0);
 693         check_domains(delta, window, graph[0]);
 694         std::cout << "finding domains' sizes\n";
 695         find_domain_sizes(graph[0], domsizes);
 696         while (check_domsizes(domsizes, domain_min_size)) {
 697             domains.push_back(a);
 698             if (DomainSearchingAlgorythm == 0) {
 699                 get_maxsized_domain(domains.back(), graph[0], domsizes);
 700             } else if (DomainSearchingAlgorythm == 1) {
 701                 get_min_domain(domains.back(), graph[0], domsizes, domain_min_size);
 702             }
 703             std::cout << "new domain: " << domains.back().size() << " atoms\n";
 704             delete_domain_from_graph(graph[0], domains.back());
 705             domsizes.resize(0);
 706             find_domain_sizes(graph[0], domsizes);
 707         }
 708         graph.erase(graph.begin());
 709         std::cout << (frnr + 1) / (window / tau_jump) - tau_jump << " window's domains have been evaluated\n";
 710         print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / (window / tau_jump) - tau_jump); // see function for details | numbers from index
 711     }
 712
 713     for (unsigned int i = 0; i < index.size(); i++) {
 714         trajectory[i] = fr.x[index[i]];
 715     }
 716     frames++;
 717
 718     #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
 719     for (unsigned int j = 0; j < bone; j++) {
 720         std::vector < RVec > temp = trajectory;
 721         MyFitNew(reference, temp, w_rls2[j], 0);
 722         for (unsigned int i = 0; i < graph.size(); i++) {
 723             update_graph(graph[i][j], temp, epsi);
 724         }
 725     }
 726     #pragma omp barrier
 727
 728     std::vector< std::vector< RVec > > tsTemp;
 729     tsTemp.resize(bone, trajectory);
 730
 731     #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
 732     for (unsigned int i = 0; i < bone; i++) {
 733         MyFitNew(reference, tsTemp[i], w_rls2[i], 0);
 734     }
 735     #pragma omp barrier
 736     testSubject.update(tsTemp, epsi, frnr);
 737     testSubject.getDomains(delta);
 738     testSubject.print(index, fnNdx_, epsi, delta, (frnr + 1) / (window / tau_jump) - tau_jump);
 739
 740     std::cout << "frame: " << frnr << " analyzed\n";
 741 }
 742
 743 void
 744 Domains::finishAnalysis(int nframes)
 745 {
 746
 747 }
 748
 749 void
 750 Domains::writeOutput()
 751 {
 752     std::cout << "\n END \n";
 753 }
 754
 755 /*! \brief
 756  * The main function for the analysis template.
 757  */
 758 int
 759 main(int argc, char *argv[])
 760 {
 761     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 762 }