src/spacetimecorr.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 <omp.h>
  46 #include <thread>
  47 #include <string>
  48 #include <algorithm>
  49
  50 #include <gromacs/trajectoryanalysis.h>
  51 #include <gromacs/pbcutil/pbc.h>
  52 #include <gromacs/pbcutil/rmpbc.h>
  53 #include <gromacs/utility/smalloc.h>
  54 #include <gromacs/math/vectypes.h>
  55 #include <gromacs/math/vec.h>
  56 #include <gromacs/math/do_fit.h>
  57
  58 #include "gromacs/analysisdata/analysisdata.h"
  59 #include "gromacs/analysisdata/modules/average.h"
  60 #include "gromacs/analysisdata/modules/histogram.h"
  61 #include "gromacs/analysisdata/modules/plot.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/options/basicoptions.h"
  64 #include "gromacs/options/filenameoption.h"
  65 #include "gromacs/options/ioptionscontainer.h"
  66 #include "gromacs/pbcutil/pbc.h"
  67 #include "gromacs/selection/selection.h"
  68 #include "gromacs/selection/selectionoption.h"
  69 #include "gromacs/trajectory/trajectoryframe.h"
  70 #include "gromacs/trajectoryanalysis/analysissettings.h"
  71 #include "gromacs/utility/arrayref.h"
  72 #include "gromacs/utility/exceptions.h"
  73 #include "gromacs/utility/stringutil.h"
  74
  75 using namespace gmx;
  76
  77 using gmx::RVec;
  78
  79 void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
  80 {
  81     FILE *file;
  82     file = std::fopen(file_name, "w+");
  83     for (int i = 1; i < trajectory.size(); i++) {
  84         std::fprintf(file, "MODEL%9d\n", i);
  85         for (int j = 0; j < trajectory[i].size(); j++) {
  86             std::fprintf(file, "ATOM  %5d   CA LYS     1    %8.3f%8.3f%8.3f  1.00  0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10);
  87         }
  88         std::fprintf(file, "ENDMDL\n");
  89     }
  90 }
  91 /*
  92 1 -  6        Record name     "ATOM  "
  93 7 - 11        Integer         Atom serial number.
  94 13 - 16        Atom            Atom name.
  95 17             Character       Alternate location indicator.
  96 18 - 20        Residue name    Residue name.
  97 22             Character       Chain identifier.
  98 23 - 26        Integer         Residue sequence number.
  99 27             AChar           Code for insertion of residues.
 100 31 - 38        Real(8.3)       Orthogonal coordinates for X in Angstroms.
 101 39 - 46        Real(8.3)       Orthogonal coordinates for Y in Angstroms.
 102 47 - 54        Real(8.3)       Orthogonal coordinates for Z in Angstroms.
 103 55 - 60        Real(6.2)       Occupancy.
 104 61 - 66        Real(6.2)       Temperature factor (Default = 0.0).
 105 73 - 76        LString(4)      Segment identifier, left-justified.
 106 77 - 78        LString(2)      Element symbol, right-justified.
 107 79 - 80        LString(2)      Charge on the atom.
 108 */
 109
 110 void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start)
 111 {
 112     FILE *file;
 113     file = std::fopen(file_name, "w+");
 114     for (int i = start; i < correlations.size(); i++) {
 115         if (correlations.size() - start > 1) {
 116             std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
 117         }
 118         if (i % 50 == 0) {
 119             std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
 120         }
 121         for (int j = 0; j < correlations[i].size(); j++) {
 122             for (int f = 0; f < correlations[i][j].size(); f++) {
 123                 std::fprintf(file, "%3.2f ", correlations[i][j][f]);
 124             }
 125             std::fprintf(file, "\n");
 126         }
 127     }
 128     std::fclose(file);
 129 }
 130
 131 void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std::vector< int > > rout, const char* file_name)
 132 {
 133     FILE *file;
 134     file = std::fopen(file_name, "w+");
 135     std::fprintf(file, "correlations >= %0.2f\n", crl_border);
 136     for (int i = 0; i < rout.size(); i++) {
 137         for (int j = 0; j < rout[i].size(); j++) {
 138             std::fprintf(file, "%3d ", indx[rout[i][j]] + 1);
 139         }
 140         std::fprintf(file, "\n");
 141     }
 142     std::fclose(file);
 143 }
 144
 145 int dive_in(std::vector< std::vector< std::pair< int, int > > > graph, int start, int steps, int frames, int depth, int depth_go, std::vector< std::pair< int, int > > path)
 146 {
 147     int steps_old = 0, steps_new = 0;
 148     if (depth > frames) {
 149         return steps;
 150     }
 151     bool flag = true;
 152     for (int i = 0; i < path.size(); i++) {
 153         if (path[i].first == start) {
 154             flag = false;
 155         }
 156     }
 157     if (flag) {
 158         path.push_back(std::make_pair(start, depth_go));
 159     }
 160     for (int i = 0; i < graph[start].size(); i++) {
 161         flag = true;
 162         for (int j = 0; j < path.size(); j++) {
 163             if (path[j].first == graph[start][i].first) {
 164                 flag = false;
 165             }
 166         }
 167         if (flag) {
 168             steps_new = dive_in(graph, graph[start][i].first, steps + 1, frames, depth + graph[start][i].second, graph[start][i].second, path);
 169             if (steps_old < steps_new) {
 170                 steps_old = steps_new;
 171             }
 172         }
 173     }
 174     return (std::max(steps_old, steps));
 175 }
 176
 177 bool mysortfunc (std::vector<int> a, std::vector<int> b) { return (a.size() > b.size()); }
 178
 179 /*! \brief
 180  * Class used to compute free volume in a simulations box.
 181  *
 182  * Inherits TrajectoryAnalysisModule and all functions from there.
 183  * Does not implement any new functionality.
 184  *
 185  * \ingroup module_trajectoryanalysis
 186  */
 187
 188
 189 class Domains : public TrajectoryAnalysisModule
 190 {
 191     public:
 192
 193         Domains();
 194         virtual ~Domains();
 195
 196         //! Set the options and setting
 197         virtual void initOptions(IOptionsContainer          *options,
 198                                  TrajectoryAnalysisSettings *settings);
 199
 200         //! First routine called by the analysis framework
 201         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 202         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 203                                   const TopologyInformation        &top);
 204
 205         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 206                                          const t_trxframe                   &fr);
 207
 208         //! Call for each frame of the trajectory
 209         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 210         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 211                                   TrajectoryAnalysisModuleData *pdata);
 212
 213         //! Last routine called by the analysis framework
 214         // virtual void finishAnalysis(t_pbc *pbc);
 215         virtual void finishAnalysis(int nframes);
 216
 217         //! Routine to write output, that is additional over the built-in
 218         virtual void writeOutput();
 219
 220     private:
 221
 222         std::string                                                 fnNdx_;
 223         SelectionList                                               sel_;
 224
 225         std::vector< std::vector< int > >                           domains;
 226         std::vector< std::vector< int > >                           domains_local;
 227         std::vector< std::vector< RVec > >                          trajectory;
 228         std::vector< std::vector< RVec > >                          frankenstein_trajectory;
 229
 230
 231         std::vector< int >                                          index;
 232         std::vector< int >                                          numbers;
 233         int                                                         frames              = 0;
 234         int                                                         basic_frame         = 0;
 235         int                                                         tau                 = 0;
 236         double                                                      crl_border          = 0;
 237         real                                                        **w_rls;
 238
 239         int                                                         domains_ePBC;
 240         // Copy and assign disallowed by base.
 241 };
 242
 243 Domains::Domains(): TrajectoryAnalysisModule()
 244 {
 245 }
 246
 247 Domains::~Domains()
 248 {
 249 }
 250
 251 void
 252 Domains::initOptions(IOptionsContainer          *options,
 253                      TrajectoryAnalysisSettings *settings)
 254 {
 255     static const char *const desc[] = {
 256         "[THISMODULE] to be done"
 257     };
 258     // Add the descriptive text (program help text) to the options
 259     settings->setHelpText(desc);
 260     // Add option for output file name
 261     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 262                             .store(&fnNdx_).defaultBasename("domains")
 263                             .description("Index file from the domains"));
 264     // Add option for tau constant
 265     options->addOption(gmx::IntegerOption("tau")
 266                             .store(&tau)
 267                             .description("number of frames for time travel"));
 268     // Add option for crl_border constant
 269     options->addOption(DoubleOption("crl")
 270                             .store(&crl_border)
 271                             .description("make graph based on corrs > constant"));
 272     // Add option for selection list
 273     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 274                            .required().dynamicMask().multiValue()
 275                            .description("Domains to form rigid skeleton"));
 276     // Control input settings
 277     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 278     settings->setPBC(true);
 279 }
 280
 281 void
 282 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 283                       const TopologyInformation        &top)
 284 {
 285     domains_ePBC = top.ePBC();
 286 }
 287
 288 void
 289 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 290                              const t_trxframe                       &fr)
 291 {
 292     ConstArrayRef< int > atomind  = sel_[0].atomIndices();
 293     index.resize(0);
 294     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 295         index.push_back(*ai);
 296     }
 297
 298     trajectory.resize(1);
 299     trajectory.back().resize(index.size());
 300
 301     for (int i = 0; i < index.size(); i++) {
 302         trajectory.back()[i] = fr.x[index[i]];
 303     }
 304
 305     domains.resize(sel_.size());
 306     for (int i = 1; i < sel_.size(); i++) {
 307         atomind  = sel_[i].atomIndices();
 308         for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 309             domains[i].push_back(*ai);
 310         }
 311     }
 312
 313     for (int i = 0; i < domains.size(); i++) {
 314         for (int j = 0; j < domains[i].size(); j++) {
 315             int k = 0;
 316             while (index[k] != domains[i][j]) {
 317                 k++;
 318             }
 319             domains[i][j] = k;
 320         }
 321     }
 322
 323     for (int i = 0; i < domains.size(); i++) {
 324         if (domains[i].size() >= 2) {
 325             domains_local.push_back(domains[i]);
 326             for (int k = 0; k < domains[i].size(); k++) {
 327                 for (int j = i + 1; j < domains.size(); j++) {
 328                     for (int x = 0; x < domains[j].size(); x++) {
 329                         if (domains[j][x] == domains[i][k]) {
 330                             domains[j].erase(domains[j].begin() + x);
 331                         }
 332                     }
 333                 }
 334             }
 335         }
 336     }
 337
 338     domains.resize(0);
 339     domains = domains_local;
 340
 341     std::vector< bool >                 flags;
 342     flags.resize(index.size(), true);
 343
 344     for (int i = 0; i < domains.size(); i++) {
 345         for (int j = 0; j < domains[i].size(); j++) {
 346             flags[domains[i][j]] = false;
 347         }
 348     }
 349
 350     for (int i = 0; i < index.size(); i++) {
 351         if (flags[i]) {
 352             int a, b, dist = 90000001;
 353             rvec temp;
 354                 for (int j = 0; j < domains.size(); j++) {
 355                     for (int k = 0; k < domains[j].size(); k++) {
 356                         rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp);
 357                         if (norm(temp) <= dist) {
 358                             dist = norm(temp);
 359                             a = j;
 360                             b = k;
 361                         }
 362                     }
 363                 }
 364             domains_local[a].push_back(i);
 365             flags[i] = false;
 366         }
 367     }
 368
 369     frankenstein_trajectory.resize(trajectory.size());
 370     frankenstein_trajectory.back() = trajectory.back();
 371
 372     frames++;
 373 }
 374
 375 void
 376 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 377                       TrajectoryAnalysisModuleData *pdata)
 378 {
 379     trajectory.resize(trajectory.size() + 1);
 380     trajectory.back().resize(index.size());
 381
 382     for (int i = 0; i < index.size(); i++) {
 383         trajectory.back()[i] = fr.x[index[i]];
 384     }
 385
 386     snew(w_rls, domains_local.size() + 1);
 387     for (int i = 0; i < domains_local.size(); i++) {
 388         snew(w_rls[i], index.size());
 389         for (int j = 0; j < index.size(); j++) {
 390             w_rls[i][j] = 0;
 391         }
 392         for (int j = 0; j < domains_local[i].size(); j++) {
 393             w_rls[i][domains_local[i][j]] = 1;
 394         }
 395     }
 396     snew(w_rls[domains_local.size()], index.size());
 397     for (int i = 0; i < index.size(); i++) {
 398         w_rls[domains_local.size()][i] = 1;
 399     }
 400
 401     frankenstein_trajectory.resize(trajectory.size());
 402     frankenstein_trajectory.back().resize(index.size());
 403
 404     for (int i = 0; i < domains_local.size(); i++) {
 405         rvec *basic, *traj;
 406         snew(basic, index.size());
 407         for (int k = 0; k < index.size(); k++) {
 408             copy_rvec(trajectory[basic_frame][k].as_vec(), basic[k]);
 409         }
 410         snew(traj, index.size());
 411         for (int k = 0; k < index.size(); k++) {
 412             copy_rvec(trajectory.back()[k].as_vec(), traj[k]);
 413         }
 414         reset_x(index.size(), NULL, index.size(), NULL, basic, w_rls[i]);
 415         reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[i]);
 416         do_fit(index.size(), w_rls[i], basic, traj);
 417
 418         for (int j = 0; j < index.size(); j++) {
 419             if (w_rls[i][j] == 0) {
 420                 copy_rvec(basic[j], traj[j]);
 421             }
 422         }
 423         reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[domains_local.size()]);
 424
 425         for (int j = 0; j < domains_local[i].size(); j++) {
 426             frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
 427         }
 428
 429         sfree(basic);
 430         sfree(traj);
 431     }
 432
 433     frames++;
 434 }
 435
 436 //spacetimecorr -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test.ndx'
 437
 438 void
 439 Domains::finishAnalysis(int nframes)
 440 {
 441     //make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/10k_frames.pdb");
 442
 443     //frankenstein_trajectory = trajectory;
 444
 445     std::vector< std::vector< std::vector< double > > > crltns;
 446     int k = 1000, m = 0;
 447     if (tau > 0) {
 448         k = tau;
 449     }
 450     crltns.resize(k + 1);
 451     for (int i = 0; i < crltns.size(); i++) {
 452         crltns[i].resize(index.size());
 453         for (int j = 0; j < crltns[i].size(); j++) {
 454             crltns[i][j].resize(index.size(), 0);
 455         }
 456     }
 457
 458     std::cout << "\n" ;
 459     #pragma omp parallel
 460     {
 461         #pragma omp for schedule(dynamic)
 462         for (int i = m; i <= k; i += 1) {
 463             std::cout << " k = " << i << " ";
 464
 465             rvec *medx, *medy, temp_zero;
 466             snew(medx, index.size());
 467             snew(medy, index.size());
 468             temp_zero[0] = 0;
 469             temp_zero[1] = 0;
 470             temp_zero[2] = 0;
 471             for (int i = 0; i < index.size(); i++) {
 472                 copy_rvec(temp_zero, medx[i]);
 473                 copy_rvec(temp_zero, medy[i]);
 474             }
 475
 476             for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
 477                 for (int f = 0; f < index.size(); f++) {
 478                     rvec temp;
 479
 480                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp);
 481                     rvec_inc(medx[f], temp);
 482                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + i][f], temp);
 483                     rvec_inc(medy[f], temp);
 484                 }
 485             }
 486
 487             for (int j = 0; j < index.size(); j++) {
 488                 rvec temp;
 489                 real temp2 = frankenstein_trajectory.size();
 490                 temp2 -= i;
 491                 temp2 = 1 / temp2;
 492                 copy_rvec(medx[j], temp);
 493                 svmul(temp2, temp, medx[j]);
 494                 copy_rvec(medy[j], temp);
 495                 svmul(temp2, temp, medy[j]);
 496             }
 497
 498             std::vector< std::vector< double > > a, b, c;
 499             a.resize(index.size());
 500             b.resize(index.size());
 501             c.resize(index.size());
 502             for (int j = 0; j < index.size(); j++) {
 503                 a[j].resize(index.size(), 0);
 504                 b[j].resize(index.size(), 0);
 505                 c[j].resize(index.size(), 0);
 506             }
 507             for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
 508                 for (int f1 = 0; f1 < index.size(); f1++) {
 509                     for (int f2 = 0; f2 < index.size(); f2++) {
 510                         rvec temp1, temp2;
 511                         rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
 512                         rvec_dec(temp1, medx[f1]);
 513                         rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + i][f2], temp2);
 514                         rvec_dec(temp2, medy[f2]);
 515                         a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
 516                         b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
 517                         c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
 518                     }
 519                 }
 520             }
 521             for (int j = 0; j < index.size(); j++) {
 522                 for (int f = 0; f < index.size(); f++) {
 523                     crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 524                 }
 525             }
 526             sfree(medx);
 527             sfree(medy);
 528         }
 529     }
 530     std::cout << "\n" ;
 531
 532     std::cout << "printing correlation's file\n" ;
 533     make_correlation_file(crltns, "corrs.txt", m);
 534     std::cout << "done\n" ;
 535
 536     for (int i1 = 0; i1 < 100; i1++) {
 537         int i5 = 0;
 538         for (int i2 = 1; i2 < crltns.size(); i2++) {
 539             for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
 540                 for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {
 541                     if (crltns[i2][i3][i4] >= (double)i1 / 100) {
 542                         i5++;
 543                     }
 544                 }
 545             }
 546         }
 547         std::cout << i1 << " - " << i5 << "   ";
 548         if (i1 % 10 == 0) {
 549             std::cout << "\n" ;
 550         }
 551     }
 552
 553     std::cout << "\n\n";
 554
 555     std::vector< std::vector< int > > graph;
 556     graph.resize(index.size());
 557     for (int i = 1; i <= k; i++) {
 558         for (int j = 0; j < index.size(); j++) {
 559             for (int f = 0; f < index.size(); f++) {
 560                 if (std::abs(crltns[i][j][f]) >= crl_border) {
 561                     bool local_flag = true;
 562                     for (int f2 = 0; f2 < graph[j].size(); f2++) {
 563                         if (graph[j][f2] == f) {
 564                             local_flag = false;
 565                         }
 566                     }
 567                     if (local_flag) {
 568                         graph[j].push_back(f);
 569                     }
 570                 }
 571             }
 572         }
 573     }
 574
 575     /*for (int i = 0; i < index.size(); i++) {
 576         if (graph[i].size() > 0) {
 577             std::cout << i << " - " << graph[i].size() << "\n";
 578         }
 579     }*/
 580
 581     std::cout << "\n";
 582     std::vector< std::vector < int > > rout_old, rout_new;
 583     bool flag = true;
 584     rout_new.resize(0);
 585     rout_old.resize(index.size());
 586     for (int i = 0; i < index.size(); i++) {
 587         rout_old[i].resize(1, i);
 588     }
 589
 590     std::cout << "building routs\n\n" ;
 591
 592     while (flag) {
 593         flag = false;
 594         std::vector < int > flag3;
 595         for (long int i = 0; i < rout_old.size(); i++) {
 596             if (graph[rout_old[i].back()].size() == 0 && rout_old[i].size() > 2) {
 597                 rout_new.push_back(rout_old[i]);
 598             } else {
 599                 bool flag222 = true;
 600                 for (long int j = 0; j < graph[rout_old[i].back()].size(); j++) {
 601                     bool flag2 = true, flag22 = true;
 602                     flag3 = rout_old[i];
 603                     for (int f = 0; f < flag3.size(); f++) {
 604                         if (flag3[f] == graph[rout_old[i].back()][j]) {
 605                             flag22 = false;
 606                         }
 607                     }
 608                     if (flag22) {
 609                         flag3.push_back(graph[rout_old[i].back()][j]);
 610                         for (int f = 0; f < rout_new.size(); f++) {
 611                             if (rout_new[f] == flag3) {
 612                                 flag2 = false;
 613                             }
 614                         }
 615                         if (flag2) {
 616                             rout_new.push_back(flag3);
 617                             flag = true;
 618                             flag222 = false;
 619
 620                         }
 621                     }
 622                 }
 623                 if (flag222 && rout_old[i].size() > 2) {
 624                     rout_new.push_back(rout_old[i]);
 625                 }
 626             }
 627         }
 628         if (rout_new.size() != 0) {
 629             std::cout << rout_old.size() << " -> " << rout_new.size() << "\n";
 630             rout_old = rout_new;
 631             rout_new.resize(0);
 632         }
 633     }
 634
 635     std::sort(rout_old.begin(), rout_old.end(), mysortfunc);
 636
 637     std::cout << "\nfinished routs\n\n" ;
 638
 639     /*graph.resize(0);
 640     graph.resize(index.size());
 641     rout_new.resize(0);
 642
 643     for (int i = 0; i < rout_old.size(); i++) {
 644         int lflag2 = false;
 645         for (int j = 0; j < rout_old[i].size() - 1; j++) {
 646             int lflag1 = true;
 647             for (int f = 0; f < graph[rout_old[i][j]].size(); f++) {
 648                 if (graph[rout_old[i][j]][f] == rout_old[i][j + 1]) {
 649                     lflag1 = false;
 650                 }
 651             }
 652             if (lflag1) {
 653                 graph[rout_old[i][j]].push_back(rout_old[i][j + 1]);
 654                 lflag2 = true;
 655             }
 656         }
 657         if (lflag2) {
 658             rout_new.push_back(rout_old[i]);
 659         }
 660     }*/
 661
 662
 663     /*std::cout << "\n test old \n";
 664     for (int i = 0; i < rout_old.size(); i++) {
 665         for (int j = 0; j < rout_old[i].size(); j++) {
 666             std::cout << index[rout_old[i][j]] << " ";
 667         }
 668         std::cout << "\n";
 669     }
 670
 671     std::cout << "\n test new \n";*/
 672
 673     rout_new.resize(0);
 674     for (int i = rout_old.size() - 1; i >= 0; i--) {
 675         bool lflag1 = true;
 676         for (int j = 0; j < i; j++) {
 677             if (rout_old[i].size() < rout_old[j].size()) {
 678                 int la = 0;
 679                 for (int f = 0; f < rout_old[j].size(); f++) {
 680                     if (rout_old[j][f] == rout_old[i][la]) {
 681                         la++;
 682                     }
 683                 }
 684                 if (la == rout_old[i].size()) {
 685                     lflag1 = false;
 686                 }
 687             }
 688         }
 689         if (lflag1) {
 690             rout_new.push_back(rout_old[i]);
 691         }
 692     }
 693
 694     make_rout_file(crl_border, index, rout_new, "routs.txt");
 695
 696     /*for (int i = 0; i < rout_new.size(); i++) {
 697         for (int j = 0; j < rout_new[i].size(); j++) {
 698             std::cout << index[rout_new[i][j]] << " ";
 699         }
 700         std::cout << "\n";
 701     }*/
 702
 703     std::cout << "\n\n\n job done \n\n\n";
 704 }
 705
 706 void
 707 Domains::writeOutput()
 708 {
 709
 710 }
 711
 712
 713 /*! \brief
 714  * The main function for the analysis template.
 715  */
 716 int
 717 main(int argc, char *argv[])
 718 {
 719     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 720 }