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