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