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
  59 #define MAX_NTRICVEC 12
  60
  61 using namespace gmx;
  62
  63 using gmx::RVec;
  64
  65 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) {
  66     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) +
  67                     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) +
  68                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  69 }
  70
  71 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) {
  72     double t01, t02, t03, t04, t05, t06, t07;
  73     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);
  74     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);
  75     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  76     t04 = sqrt(t01 + t02 + t03);
  77     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);
  78     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  79     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);
  80     F += t04;
  81     Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  82     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);
  83     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);
  84     Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  85             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 +
  86             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);
  87     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 +
  88             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 -
  89             2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  90     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 -
  91             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);
  92 }
  93
  94 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
  95     double t0 = 0, t1 = 0, t2 = 0;
  96     for (unsigned int i = 0; i < b.size(); i++) {
  97         t0 = static_cast< double >(b[i][0]);
  98         t1 = static_cast< double >(b[i][1]);
  99         t2 = static_cast< double >(b[i][2]);
 100         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));
 101         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));
 102         b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
 103     }
 104 }
 105
 106 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
 107     midA[0] = 0;
 108     midA[1] = 0;
 109     midA[2] = 0;
 110
 111     midB[0] = 0;
 112     midB[1] = 0;
 113     midB[2] = 0;
 114
 115     for (unsigned int i = 0; i < pairs.size(); i++) {
 116         rvec_inc(midA, a[pairs[i].first]);
 117         rvec_inc(midB, b[pairs[i].second]);
 118     }
 119     midA[0] /= pairs.size();
 120     midA[1] /= pairs.size();
 121     midA[2] /= pairs.size();
 122
 123     midB[0] /= pairs.size();
 124     midB[1] /= pairs.size();
 125     midB[2] /= pairs.size();
 126 }
 127
 128 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
 129     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 130     RVec ma, mb;
 131     CalcMid(a, b, ma, mb, pairs);
 132     rvec_dec(ma, mb);
 133     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;
 134     for (unsigned int i = 0; i < pairs.size(); i++) {
 135         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]),
 136                 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);
 137     }
 138     if (FtCnst == 0) {
 139         FtCnst = 0.00001;
 140     }
 141     if (f1 == 0) {
 142         return;
 143     } else {
 144         int count = 0;
 145         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
 146             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
 147             for (unsigned int i = 0; i < pairs.size(); i++) {
 148                 searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
 149                                static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 150                                static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
 151                                static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 152             }
 153             while (f2 != f1) {
 154                 f2 = 0;
 155                 for (unsigned int i = 0; i < pairs.size(); i++) {
 156                     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]),
 157                             static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
 158                             static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
 159                 }
 160                 count++;
 161                 if (f2 < f1) {
 162                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 163                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 164                     break;
 165                 } else {
 166                     //l *= 0.85;
 167                     l *= 0.50;
 168                     /*if (count >= 14) {
 169                         std::cout << count << " " << l << "\n";
 170                     }*/
 171                     if (l == DBL_MIN) { /*DBL_TRUE_MIN*/
 172                         break;
 173                     }
 174                 }
 175             }
 176             count++;
 177             if (count % 100000 == 0) {
 178                 std::cout << "+100k\n";
 179             }
 180         }
 181         ApplyFit(b, x, y, z, A, B, C);
 182     }
 183 }
 184
 185 struct node {
 186     short int n;
 187     RVec r;
 188     bool yep;
 189 };
 190
 191 void make_graph(unsigned long mgwi_natoms, std::vector < RVec > mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
 192 {
 193     mgwi_graph.resize(mgwi_natoms);
 194     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 195         mgwi_graph[i].resize(mgwi_natoms);
 196     }
 197     for (unsigned int i = 0; i < mgwi_natoms; i++) {
 198         for (unsigned int j = 0; j < mgwi_natoms; j++) {
 199             rvec_sub(mgwi_x[i].as_vec(), mgwi_x[j].as_vec(), mgwi_graph[i][j].r.as_vec());
 200             mgwi_graph[i][j].n = 0;
 201         }
 202     }
 203 }
 204
 205 void update_graph(std::vector< std::vector< node > > &ugwi_graph, std::vector < RVec > ugwi_x, /*long*/ double ugwi_epsi) {
 206     RVec ugwi_temp;
 207     for (unsigned int i = 0; i < ugwi_graph.size(); i++) {
 208         for (unsigned int j = i; j < ugwi_graph.size(); j++) {
 209             rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), ugwi_temp.as_vec());
 210             rvec_dec(ugwi_temp.as_vec(), ugwi_graph[i][j].r.as_vec());
 211             if (static_cast< double >(norm(ugwi_temp)) <= ugwi_epsi) {
 212                 if (i == j) {
 213                     ugwi_graph[i][j].n++;
 214                 }
 215                 else {
 216                     ugwi_graph[i][j].n++;
 217                     ugwi_graph[j][i].n++;
 218                 }
 219             }
 220         }
 221     }
 222 }
 223
 224 void check_domains(double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
 225     for (unsigned int k = 0; k < cd_graph.size(); k++) {
 226         for (unsigned int i = 0; i < cd_graph[1].size(); i++) {
 227             for (unsigned int j = 0; j < cd_graph[1].size(); j++) {
 228                 if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
 229                     cd_graph[k][i][j].yep = true;
 230                 }
 231                 else {
 232                     cd_graph[k][i][j].yep = false;
 233                 }
 234             }
 235         }
 236     }
 237 }
 238
 239 void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
 240     fds_domsizes.resize(fds_graph.size());
 241     for (unsigned int i = 0; i < fds_graph.size(); i++) {
 242         fds_domsizes[i].resize(fds_graph[1].size(), 0);
 243         for (unsigned int j = 0; j < fds_graph[1].size(); j++) {
 244             for (unsigned int k = 0; k < fds_graph[1].size(); k++) {
 245                 if (fds_graph[i][j][k].yep) {
 246                     fds_domsizes[i][j]++;
 247                 }
 248             }
 249         }
 250     }
 251 }
 252
 253 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) {
 254     unsigned int gmd_number1 = 0, gmd_number2 = 0;
 255     for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
 256         for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
 257             if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
 258                 gmd_number1 = i;
 259                 gmd_number2 = j;
 260             }
 261         }
 262     }
 263     gmd_max_d.resize(0);
 264     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 265         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 266             gmd_max_d.push_back(i);
 267         }
 268     }
 269 }
 270
 271 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) {
 272     unsigned int gmd_number1 = 0, gmd_number2 = 0;
 273     for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
 274         for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
 275             if (gmd_domsizes[gmd_number1][gmd_number2] < gmd_min_dom_size && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 276                 gmd_number1 = i;
 277                 gmd_number2 = j;
 278             }
 279             if (gmd_domsizes[i][j] <= gmd_domsizes[gmd_number1][gmd_number2] && gmd_domsizes[i][j] >= gmd_min_dom_size) {
 280                 gmd_number1 = i;
 281                 gmd_number2 = j;
 282             }
 283         }
 284     }
 285     gmd_min_d.resize(0);
 286     for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
 287         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 288             gmd_min_d.push_back(i);
 289         }
 290     }
 291 }
 292
 293 void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< unsigned int > ddf_domain) {
 294     for (unsigned int i = 0; i < ddf_domain.size(); i++) {
 295         for (unsigned int k = 0; k < ddf_graph.size(); k++) {
 296             for (unsigned int j = 0; j < ddf_graph[1].size(); j++) {
 297                 if (ddf_graph[k][ddf_domain[i]][j].yep) {
 298                     ddf_graph[k][ddf_domain[i]][j].yep = false;
 299                 }
 300                 if (ddf_graph[k][j][ddf_domain[i]].yep) {
 301                     ddf_graph[k][j][ddf_domain[i]].yep = false;
 302                 }
 303             }
 304         }
 305     }
 306 }
 307
 308 bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
 309     for (unsigned int i = 0; i < cd_domsizes.size(); i++) {
 310         for (unsigned int j = 0; j < cd_domsizes[0].size(); j++) {
 311             if (cd_domsizes[i][j] >= cd_domain_min_size) {
 312                 return true;
 313             }
 314         }
 315     }
 316     return false;
 317 }
 318
 319 void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta) {
 320     FILE                *fpNdx_;
 321     fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
 322     int write_count;
 323     for (unsigned int i = 0; i < pd_domains.size(); i++) {
 324         std::fprintf(fpNdx_, "[domain_№_%3d_%3d_%4.3f_%4.3f]\n", i + 1, dms, epsi, delta);
 325         write_count = 0;
 326         for (unsigned int j = 0; j < pd_domains[i].size(); j++) {
 327             write_count++;
 328             if (write_count > 20) {
 329                 write_count -= 20;
 330                 std::fprintf(fpNdx_, "\n");
 331             }
 332             std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
 333             std::printf("%5d ", index[pd_domains[i][j]] + 1);
 334         }
 335         std::printf("\n\n");
 336         std::fprintf(fpNdx_,"\n\n");
 337     }
 338     std::fprintf(fpNdx_,"\n");
 339     std::fclose(fpNdx_);
 340 }
 341
 342 /*! \brief
 343  * \ingroup module_trajectoryanalysis
 344  */
 345 class Domains : public TrajectoryAnalysisModule
 346 {
 347     public:
 348
 349         Domains();
 350         virtual ~Domains();
 351
 352         //! Set the options and setting
 353         virtual void initOptions(IOptionsContainer          *options,
 354                                  TrajectoryAnalysisSettings *settings);
 355
 356         //! First routine called by the analysis framework
 357         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 358         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 359                                   const TopologyInformation        &top);
 360
 361         //! Call for each frame of the trajectory
 362         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 363         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 364                                   TrajectoryAnalysisModuleData *pdata);
 365
 366         //! Last routine called by the analysis framework
 367         // virtual void finishAnalysis(t_pbc *pbc);
 368         virtual void finishAnalysis(int nframes);
 369
 370         //! Routine to write output, that is additional over the built-in
 371         virtual void writeOutput();
 372
 373     private:
 374
 375         std::string                                                 fnNdx_;
 376
 377         //std::vector< std::vector< std::vector< node > > >           graph;
 378         std::vector< std::vector< std::vector< std::vector< node > > > >    graph;
 379
 380
 381         std::vector< std::vector< unsigned int > >                  domains;
 382         std::vector< std::vector< int > >                           domsizes;
 383
 384         std::vector< int >                                          index;
 385         std::vector< int >                                          numbers;
 386         std::vector< RVec >                                         trajectory;
 387         std::vector< RVec >                                         reference;
 388         Selection                                                   selec;
 389         int                                                         frames              = 0;
 390         int                                                         domain_min_size     = 4; // selectable
 391
 392         std::vector< std::vector< std::pair< unsigned int, unsigned int > > >         w_rls2;
 393         unsigned int                                                bone;
 394         double                                                      delta               = 0.95; // selectable
 395         double                                                      epsi                = 0.10; // selectable
 396
 397         int                                                         DomainSearchingAlgorythm = 0; // selectable
 398         // Copy and assign disallowed by base.
 399 };
 400
 401 Domains::Domains(): TrajectoryAnalysisModule()
 402 {
 403 }
 404
 405 Domains::~Domains()
 406 {
 407 }
 408
 409 void
 410 Domains::initOptions(IOptionsContainer          *options,
 411                      TrajectoryAnalysisSettings *settings)
 412 {
 413     static const char *const desc[] = {
 414         "[THISMODULE] to be done"
 415     };
 416     // Add the descriptive text (program help text) to the options
 417     settings->setHelpText(desc);
 418     // Add option for selecting a subset of atoms
 419     options->addOption(SelectionOption("select")
 420                            .store(&selec).required()
 421                            .description("Atoms that are considered as part of the excluded volume"));
 422     // Add option for output file name
 423     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 424                             .store(&fnNdx_).defaultBasename("domains")
 425                             .description("Index file from the domains"));
 426     // Add option for domain min size constant
 427     options->addOption(gmx::IntegerOption("dms")
 428                             .store(&domain_min_size)
 429                             .description("minimum domain size, should be >= 4"));
 430     // Add option for Domain's Searching Algorythm
 431     options->addOption(gmx::IntegerOption("DSA")
 432                             .store(&DomainSearchingAlgorythm)
 433                             .description("Domain's Searching Algorythm: 0 == default (from bigger to smaller) | 1 == (from smaller to bigger)"));
 434     // Add option for epsi constant
 435     options->addOption(DoubleOption("epsilon")
 436                             .store(&epsi)
 437                             .description("thermal vibrations' constant"));
 438     // Add option for delta constant
 439     options->addOption(DoubleOption("delta")
 440                             .store(&delta)
 441                             .description("domain membership probability"));
 442     // Control input settings
 443     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 444     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 445     settings->setPBC(true);
 446 }
 447
 448 void
 449 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 450                       const TopologyInformation        &top)
 451 {
 452     index.resize(0);
 453     for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
 454         index.push_back(*ai);
 455     }
 456
 457     reference.resize(0);
 458     if (top.hasFullTopology()) {
 459         for (unsigned int i = 0; i < index.size(); i++) {
 460             reference.push_back(top.x().at(index[i]));
 461         }
 462     }
 463
 464     graph.resize(0);
 465
 466     bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
 467     graph[0].resize(bone);
 468     w_rls2.resize(bone);
 469     for (unsigned int i = 0; i < bone; i++) {
 470         w_rls2[i].resize(0);
 471         for (unsigned int j = 0; j < static_cast< unsigned int >(domain_min_size); j++) {
 472             w_rls2[i].push_back(std::make_pair(j + i, j + i));
 473         }
 474         make_graph(index.size(), reference, graph[0][i]);
 475     }
 476
 477     trajectory.resize(index.size());
 478 }
 479
 480 void
 481 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 482 {
 483     for (unsigned int i = 0; i < index.size(); i++) {
 484         trajectory[i] = fr.x[index[i]];
 485     }
 486     frames++;
 487
 488     #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi)
 489     for (unsigned int j = 0; j < bone; j++) {
 490         std::vector < RVec > temp = trajectory;
 491         MyFitNew(reference, temp, w_rls2[j], 0);
 492         update_graph(graph[0][j], temp, epsi);
 493     }
 494     #pragma omp barrier
 495     std::cout << "frame: " << frames << " analyzed\n";
 496 }
 497
 498 void
 499 Domains::finishAnalysis(int nframes)
 500 {
 501     frames -= 1;
 502
 503     std::cout << "final cheking\n";
 504     check_domains(delta, frames, graph[0]);
 505
 506     std::cout << "finding domains' sizes\n";
 507     find_domain_sizes(graph[0], domsizes);
 508
 509     std::cout << "finding domains\n";
 510     std::vector< unsigned int > a;
 511     a.resize(0);
 512     while (check_domsizes(domsizes, domain_min_size)) {
 513         domains.push_back(a);
 514         if (DomainSearchingAlgorythm == 0) {
 515             get_maxsized_domain(domains.back(), graph[0], domsizes);
 516         } else if (DomainSearchingAlgorythm == 1) {
 517             get_min_domain(domains.back(), graph[0], domsizes, domain_min_size);
 518         }
 519         std::cout << "new domain: " << domains.back().size() << " atoms\n";
 520         delete_domain_from_graph(graph[0], domains.back());
 521         domsizes.resize(0);
 522         find_domain_sizes(graph[0], domsizes);
 523     }
 524 }
 525
 526 void
 527 Domains::writeOutput()
 528 {
 529     std::cout << "making output file\n";
 530     print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta); // see function for details | numbers from index
 531     std::cout << "\n END \n";
 532 }
 533
 534 /*! \brief
 535  * The main function for the analysis template.
 536  */
 537 int
 538 main(int argc, char *argv[])
 539 {
 540     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 541 }