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