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
  44 #include <iostream>
  45 #include <algorithm>
  46 #include <vector>
  47 #include <math.h>
  48 #include <cmath>
  49 #include <omp.h>
  50 #include <string>
  51 #include <sstream>
  52
  53 #include <gromacs/trajectoryanalysis.h>
  54 #include <gromacs/utility/smalloc.h>
  55 #include <gromacs/math/do_fit.h>
  56 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  57
  58 using namespace gmx;
  59
  60 using gmx::RVec;
  61
  62 const std::vector< std::vector< std::vector< double > > > read_correlation_matrix_file(const char* file_name, unsigned long size, unsigned int length)
  63 {
  64     FILE *file;
  65     std::vector< std::vector< std::vector< double > > > crrlts;
  66     file = std::fopen(file_name, "r+");
  67     std::vector< double > c;
  68     c.resize(size, 0);
  69     std::vector< std::vector< double > > b;
  70     b.resize(size, c);
  71     crrlts.resize(length + 1, b);
  72     char a[100];
  73     std::cout << "reading correlations from a file:\n";
  74     for (unsigned int i = 0; i < length + 1; i++) {
  75         //char *t1 = std::fgets(a, 100, file);
  76         int t0 = 0, t1 = std::fscanf(file, "%d\n", &t0);
  77         std::cout << t0 << " | ";
  78         if (i % 100 == 0) {
  79             std::cout << "\n";
  80         }
  81         for (unsigned int j = 0; j < size; j++) {
  82             for (unsigned int f = 0; f < size; f++) {
  83                 int t2 = std::fscanf(file, "%lf ", &crrlts[i][j][f]);
  84             }
  85         }
  86     }
  87     std::fclose(file);
  88     std::cout << "\n";
  89     return crrlts;
  90 }
  91
  92 void make_correlation_matrix_file(const std::vector< std::vector< std::vector< double > > > correlations, const char* file_name)
  93 {
  94     FILE *file;
  95     file = std::fopen(file_name, "w+");
  96     std::cout << "writing correlation matrixes: \n";
  97     for (unsigned int i = 0; i < correlations.size(); i++) {
  98         std::fprintf(file, "%d\n", i);
  99         for (unsigned int j = 0; j < correlations[i].size(); j++) {
 100             for (unsigned int f = 0; f < correlations[i][j].size(); f++) {
 101                 std::fprintf(file, "%.6f ", correlations[i][j][f]); //~16
 102             }
 103             std::fprintf(file, "\n");
 104         }
 105         std::cout << i << " | ";
 106     }
 107     std::fclose(file);
 108     std::cout << "\n";
 109 }
 110
 111 void make_correlation_pairs_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name)
 112 {
 113     FILE *file;
 114     file = std::fopen(file_name, "w+");
 115     for (unsigned int i = 0; i < correlations.front().size(); i++) {
 116         for (unsigned int j = 0; j < correlations.front().front().size(); j++) {
 117             //std::fprintf(file, "correlation between point '%d' and point '%d'\n", i, j);
 118             std::fprintf(file, "%d %d\n", i, j);
 119             for (unsigned int k = 0; k < correlations.size(); k++) {
 120                 std::fprintf(file, "%3.4f ", correlations[k][i][j]);
 121             }
 122             std::fprintf(file, "\n");
 123         }
 124         if (i % 10 == 0) {
 125             std::cout << "correlations in row " << i << " out of " << correlations.front().size() << " completed\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< std::pair< unsigned int, unsigned 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\n", crl_border);
 136     for (unsigned int i = 0; i < rout.size(); i++) {
 137         for (unsigned int j = 0; j < rout[i].size(); j++) {
 138             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*/);
 139         }
 140         std::fprintf(file, "\n\n");
 141     }
 142     std::fclose(file);
 143 }
 144
 145 void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations,
 146                               std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout_pairs,
 147                               std::vector< int > indx,
 148                               const char* file_name)
 149 {
 150     FILE *file;
 151     file = std::fopen(file_name, "w+");
 152     for (unsigned int i = 0; i < rout_pairs.size(); i++) {
 153         for (unsigned int j = 0; j < rout_pairs[i].size(); j++) {
 154             std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first]/* + 1*/, indx[rout_pairs[i][j].second]/* + 1*/);
 155             for (unsigned int k = 0; k < correlations.size(); k++) {
 156                 std::fprintf(file, "%3.5f ", correlations[k][rout_pairs[i][j].first][rout_pairs[i][j].second]);
 157             }
 158             std::fprintf(file, "\n");
 159         }
 160     }
 161     std::fclose(file);
 162 }
 163
 164 void make_diffusion_file(const char* file_name, std::vector< double > D)
 165 {
 166     FILE *file;
 167     file = std::fopen(file_name, "w+");
 168     for (unsigned  int i = 0; i < D.size(); i++) {
 169         std::fprintf(file, "%f ", D[i]);
 170     }
 171     std::fclose(file);
 172 }
 173
 174 bool mysortfunc (std::vector< int > a, std::vector< int > b) {
 175     return (a.size() > b.size());
 176 }
 177
 178 bool isitsubset (std::vector< int > a, std::vector< int > b) {
 179     if (b.size() == 0) {
 180         return true;
 181     }
 182     std::sort(a.begin(), a.end());
 183     std::sort(b.begin(), b.end());
 184     unsigned int k = 0;
 185     for (unsigned int i = 0; i < a.size(); i++) {
 186         if (b[k] == a[i]) {
 187             k++;
 188         }
 189     }
 190     if (k == b.size()) {
 191         return true;
 192     } else {
 193         return false;
 194     }
 195 }
 196
 197 const std::vector< std::vector< std::vector< double > > > correlation_evaluation(std::vector< RVec > ref, std::vector< std::vector< RVec > > traj, unsigned int tauE) {
 198     std::vector< std::vector< std::vector< double > > > crl;
 199     crl.resize(tauE + 1);
 200     for (unsigned int i = 0; i < crl.size(); i++) {
 201         crl[i].resize(traj.front().size());
 202         for (unsigned int j = 0; j < crl[i].size(); j++) {
 203             crl[i][j].resize(traj.front().size(), 0);
 204         }
 205     }
 206     RVec temp_zero;
 207     temp_zero[0] = 0;
 208     temp_zero[1] = 0;
 209     temp_zero[2] = 0;
 210
 211     std::vector< double > d;
 212     d.resize(traj.front().size(), 0);
 213
 214     #pragma omp parallel for schedule(dynamic)
 215     for (unsigned int i = 0; i <= tauE; i++) {
 216         RVec temp1, temp2;
 217         std::vector< std::vector< double > > a, b, c;
 218         a.resize(traj.front().size(), d);
 219         b.resize(traj.front().size(), d);
 220         c.resize(traj.front().size(), d);
 221         for (unsigned int j = 0; j < traj.size() - i - 1; j++) {
 222             for (unsigned int f1 = 0; f1 < traj.front().size(); f1++) {
 223                 for (unsigned int f2 = 0; f2 < traj.front().size(); f2++) {
 224                     temp1 = traj[j][f1] - ref[f1];
 225                     temp2 = traj[j + i][f2] - ref[f2];
 226                     a[f1][f2] +=    (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
 227                                     static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
 228                                     static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
 229                     b[f1][f2] +=    (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
 230                                     static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
 231                                     static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
 232                     c[f1][f2] +=    (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
 233                                     static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
 234                                     static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
 235                 }
 236             }
 237         }
 238         for (unsigned int j = 0; j < traj.front().size(); j++) {
 239             for (unsigned int f = 0; f < traj.front().size(); f++) {
 240                 crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 241             }
 242         }
 243         std::cout << i << " | ";
 244         if (i % 100 == 0) {
 245             std::cout << "\n";
 246         }
 247     }
 248     #pragma omp barrier
 249     std::cout << "\n";
 250     return crl;
 251 }
 252
 253 void graph_calculation( std::vector< std::vector< std::pair< double, long int > > > &graph, std::vector< std::vector< unsigned int > > &s_graph,
 254                         std::vector< std::vector< std::pair< unsigned int, unsigned int > > > &s_graph_rbr,
 255                         std::vector< std::vector< RVec > > traj, std::vector< RVec > ref,
 256                         const std::vector< std::vector< std::vector<double > > > crl, double crl_b, double e_rad, unsigned int tauB, unsigned int tauE) {
 257     graph.resize(traj.front().size());
 258     for (unsigned int i = 0; i < traj.front().size(); i++) {
 259         graph[i].resize(traj.front().size(), std::make_pair(0, -1));
 260     }
 261     RVec temp;
 262     for (unsigned int i = tauB; i <= tauE; i++) {
 263         for (unsigned int j = 0; j < traj.front().size(); j++) {
 264             for (unsigned int f = j; f < traj.front().size(); f++) {
 265                 temp = ref[i] - ref[f];
 266                 if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b &&
 267                     static_cast<double>(norm(temp)) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
 268                     if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
 269                         graph[j][f].first = crl[i][j][f];
 270                     } else {
 271                         graph[j][f].first = crl[i][f][j];
 272                     }
 273                     graph[j][f].second = i;
 274                 }
 275             }
 276         }
 277     }
 278     std::cout << "crl analysed\n";
 279     std::vector< bool > graph_flags;
 280     graph_flags.resize(traj.front().size(), true);
 281     std::vector< unsigned int > a;
 282     a.resize(0);
 283     std::vector< std::pair< unsigned int, unsigned int > > b;
 284     b.resize(0);
 285     std::vector< unsigned int > que1, que2, que3;
 286     for (unsigned int i = 0; i < traj.front().size(); i++) {
 287         if (graph_flags[i]) {
 288             s_graph.push_back(a);
 289             s_graph_rbr.push_back(b);
 290             que1.resize(0);
 291             que2.resize(0);
 292             que3.resize(0);
 293             que1.push_back(i);
 294             que3.push_back(i);
 295             graph_flags[i] = false;
 296             while(que1.size() > 0) {
 297                 que2.clear();
 298                 for (unsigned int k = 0; k < que1.size(); k++) {
 299                     for (unsigned int j = 0; j < traj.front().size(); j++) {
 300                         if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
 301                             que2.push_back(j);
 302                             graph_flags[j] = false;
 303                         }
 304                     }
 305                 }
 306                 que1 = que2;
 307                 for (unsigned int j = 0; j < que2.size(); j++) {
 308                     que3.push_back(que2[j]);
 309                 }
 310             }
 311             s_graph.back() = que3;
 312             for (unsigned int j = 0; j < que3.size(); j++) {
 313                 for (unsigned int k = 0; k < traj.front().size(); k++) {
 314                     if (graph[que3[j]][k].second > -1) {
 315                         s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
 316                     }
 317                 }
 318             }
 319             //std::cout << s_graph.back().size() << "   ";
 320         }
 321     }
 322 }
 323
 324 bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
 325     return i.second < j.second;
 326 }
 327
 328 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< unsigned int, unsigned int > > > &rout_n, unsigned long indxSize,
 329                                 std::vector< std::vector< std::pair< double, long > > > graph, std::vector< std::vector< unsigned int > > s_graph,
 330                                 std::vector< std::vector< std::pair< unsigned int, unsigned int > > > s_graph_rbr) {
 331     std::vector< double > key;
 332     std::vector< long > path;
 333     std::vector< std::pair< unsigned int, double > > que;
 334     std::vector< std::pair< unsigned int, unsigned int > > a;
 335     unsigned int v;
 336     for (unsigned int i = 0; i < s_graph.size(); i++) {
 337         key.resize(0);
 338         path.resize(0);
 339         que.resize(0);
 340         v = 0;
 341         if (s_graph[i].size() > 2) {
 342             key.resize(indxSize, 2);
 343             path.resize(indxSize, -1);
 344             key[s_graph[i][0]] = 0;
 345             for (unsigned int j = 0; j < s_graph[i].size(); j++) {
 346                 que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
 347             }
 348             std::sort(que.begin(), que.end(), myfunction);
 349             while (!que.empty()) {
 350                 v = que[0].first;
 351                 que.erase(que.begin());
 352                 for (unsigned int j = 0; j < s_graph_rbr[i].size(); j++) {
 353                     long u = -1;
 354                     if (s_graph_rbr[i][j].first == v) {
 355                         u = s_graph_rbr[i][j].second;
 356                     } else if (s_graph_rbr[i][j].second == v) {
 357                         u = s_graph_rbr[i][j].first;
 358                     }
 359                     bool flag = false;
 360                     unsigned long pos = 0;
 361                     for (unsigned long k = 0; k < que.size(); k++) {
 362                         if (que[k].first == u) {
 363                             flag = true;
 364                             pos = k;
 365                             k = que.size() + 1;
 366                         }
 367                     }
 368                     if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
 369                         path[static_cast< unsigned long >(u)] = v;
 370                         key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
 371                         que[pos].second = key[static_cast< unsigned long >(u)];
 372                         sort(que.begin(), que.end(), myfunction);
 373                     }
 374                 }
 375             }
 376             a.resize(0);
 377             rout_n.push_back(a);
 378             for (unsigned int j = 0; j < indxSize; j++) {
 379                 if (path[j] != -1) {
 380                     rout_n.back().push_back(std::make_pair(j, path[j]));
 381                 }
 382             }
 383         }
 384     }
 385 }
 386
 387 gmx::RVec evaluate_com(std::vector< RVec > frame) {
 388     RVec temp;
 389     temp[0] = 0;
 390     temp[1] = 0;
 391     temp[2] = 0;
 392     for (unsigned int i = 0; i < frame.size(); i++) {
 393         temp += frame[i];
 394     }
 395     temp[0] /= frame.size();
 396     temp[1] /= frame.size();
 397     temp[2] /= frame.size();
 398     return temp;
 399 }
 400
 401 void evaluate_diffusion(std::vector< std::vector< RVec > > trj, std::vector< double > &D/*, int max_frame_depth*/) {
 402     D.resize(static_cast< unsigned int >(trj.size() * 0.9 - 1)/*max_frame_depth*/);
 403     double temp;
 404     for (unsigned int i = 1; i < trj.size() * 0.9 /*max_frame_depth*/; i++) {
 405         temp = 0;
 406         for (unsigned int j = 0; j < trj.size() - 1 - i /*trj.size() - 1 - max_frame_depth*/; j++) {
 407             temp += static_cast< unsigned int >((evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2());
 408             //D[i][j] = (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2() / (2 * i);
 409         }
 410         D[i - 1] = temp / (trj.size() - 1 - i) / (2 * i * 0.000001);
 411     }
 412 }
 413
 414 /*! \brief
 415  * \ingroup module_trajectoryanalysis
 416  */
 417
 418 class SpaceTimeCorr : public TrajectoryAnalysisModule
 419 {
 420     public:
 421
 422         SpaceTimeCorr();
 423         virtual ~SpaceTimeCorr();
 424
 425         //! Set the options and setting
 426         virtual void initOptions(IOptionsContainer          *options,
 427                                  TrajectoryAnalysisSettings *settings);
 428
 429         //! First routine called by the analysis framework
 430         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 431         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 432                                   const TopologyInformation        &top);
 433
 434         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 435                                          const t_trxframe                   &fr);
 436
 437         //! Call for each frame of the trajectory
 438         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 439         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 440                                   TrajectoryAnalysisModuleData *pdata);
 441
 442         //! Last routine called by the analysis framework
 443         // virtual void finishAnalysis(t_pbc *pbc);
 444         virtual void finishAnalysis(int nframes);
 445
 446         //! Routine to write output, that is additional over the built-in
 447         virtual void writeOutput();
 448
 449     private:
 450
 451         SelectionList                                               sel_;
 452
 453         std::vector< std::vector< RVec > >                          trajectory;
 454         std::vector< RVec >                                         reference;
 455
 456         std::vector< int >                                          index;
 457         int                                                         frames              = 0;
 458         int                                                         tau                 = 0;    // selectable
 459         float                                                       crl_border          = 0;    // selectable
 460         float                                                       eff_rad             = 1.5;  // selectable
 461         std::string                                                 OutPutName;                 // selectable
 462         int                                                         mode                = 0;    // selectable
 463         std::string                                                 MtrxNm;                     // selectable
 464         // Copy and assign disallowed by base.
 465 };
 466
 467 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
 468 {
 469 }
 470
 471 SpaceTimeCorr::~SpaceTimeCorr()
 472 {
 473 }
 474
 475 void
 476 SpaceTimeCorr::initOptions(IOptionsContainer          *options,
 477                      TrajectoryAnalysisSettings *settings)
 478 {
 479     static const char *const desc[] = {
 480         "[THISMODULE] to be done"
 481     };
 482     // Add the descriptive text (program help text) to the options
 483     settings->setHelpText(desc);
 484     // Add option for output file name
 485     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 486     //                        .store(&fnNdx_).defaultBasename("domains")
 487     //                        .description("Index file from the domains"));
 488     // Add option for working mode
 489     options->addOption(gmx::IntegerOption("mode")
 490                             .store(&mode)
 491                             .description("default 0 | rdy correlation matrixes 1, need extra params"));
 492     // Add option for Matrix Input file names
 493     options->addOption(StringOption("Mtrx_in_put")
 494                             .store(&MtrxNm)
 495                             .description("mandatory if work mode == 1"));
 496     // Add option for tau constant
 497     options->addOption(gmx::IntegerOption("tau")
 498                             .store(&tau)
 499                             .description("number of frames for time travel"));
 500     // Add option for crl_border constant
 501     options->addOption(FloatOption("crl")
 502                             .store(&crl_border)
 503                             .description("make graph based on corrs > constant"));
 504     // Add option for eff_rad constant
 505     options->addOption(FloatOption("ef_rad")
 506                             .store(&eff_rad)
 507                             .description("effective radius for atoms to evaluate corrs"));
 508     // Add option for selection list
 509     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 510                            .required().dynamicMask().multiValue()
 511                            .description("Domains to form rigid skeleton"));
 512     // Add option for output file names
 513     options->addOption(StringOption("out_put")
 514                             .store(&OutPutName)
 515                             .description("<your name here> + <local file tag>.txt"));
 516     // Control input settings
 517     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 518     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 519     settings->setPBC(true);
 520 }
 521
 522 void
 523 SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 524 {
 525     ArrayRef< const int > atomind = sel_[0].atomIndices();
 526     index.resize(0);
 527     for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 528         index.push_back(*ai);
 529     }
 530     trajectory.resize(0);
 531
 532     reference.resize(0);
 533     if (top.hasFullTopology()) {
 534         for (unsigned int i = 0; i < index.size(); i++) {
 535             reference.push_back(top.x().at(index[i]));
 536         }
 537     }
 538 }
 539
 540 void
 541 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
 542 {
 543 }
 544
 545 // -s '*.pdb' -f '*.pdb' -n '*.ndx' -sf 'name' -tau 1000 -crl 0.20 -ef_rad 9000 -out_put OLA
 546 // -s '*.tpr' -f '*.xtc' -n '*.ndx' -sf 'name' -tau 1000 -crl 0.30 -ef_rad 20 -out_put 'test_run'
 547
 548 void
 549 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 550 {
 551     trajectory.resize(trajectory.size() + 1);
 552     trajectory.back().resize(index.size());
 553     for (unsigned int i = 0; i < index.size(); i++) {
 554         trajectory.back()[i] = fr.x[index[i]];
 555     }
 556     frames++;
 557 }
 558
 559 void
 560 SpaceTimeCorr::finishAnalysis(int nframes)
 561 {
 562     std::vector< std::vector< std::vector< double > > > crltns, crltns_test;
 563     std::vector< std::vector< std::pair< double, long int > > > graph;
 564     std::vector< std::vector< unsigned int > > sub_graph;
 565     std::vector< std::vector< std::pair< unsigned int, unsigned int > > > sub_graph_rbr;
 566     std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout_new;
 567     unsigned int k = 1000, m = 0;
 568     if (tau > 0) {
 569         k = static_cast< unsigned int>(tau);
 570     }
 571
 572     if (mode == 0) {
 573         std::cout << "\nCorrelation's evaluation - start\n";
 574
 575         crltns = correlation_evaluation(reference, trajectory, k);
 576
 577         for (unsigned int i1 = 0; i1 < crltns.size(); i1++) {
 578             for (unsigned int i2 = 0; i2 < crltns[i1].size(); i2++) {
 579                 for (unsigned int i3 = 0; i3 < crltns[i1][i2].size(); i3++) {
 580                     crltns[i1][i2][i3] = std::round(crltns[i1][i2][i3] * 10000) / 10000;
 581                 }
 582             }
 583         }
 584
 585         make_correlation_matrix_file(crltns, (OutPutName + "_matrix.txt").c_str());
 586         std::cout << "corelation matrix printed\n";
 587
 588         //make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
 589         //std::cout << "corelation pairs printed\n";
 590
 591         std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
 592     } else if (mode == 1) {
 593         crltns = read_correlation_matrix_file((MtrxNm).c_str(), index.size(), k);
 594
 595         for (unsigned int i1 = 0; i1 < crltns.size(); i1++) {
 596             for (unsigned int i2 = 0; i2 < crltns[i1].size(); i2++) {
 597                 for (unsigned int i3 = 0; i3 < crltns[i1][i2].size(); i3++) {
 598                     crltns[i1][i2][i3] = std::round(crltns[i1][i2][i3] * 10000) / 10000;
 599                 }
 600             }
 601         }
 602     }
 603
 604
 605     /*
 606     int ola001 = 0;
 607     for (int i = 0; i < crltns.size(); i++) {
 608         for (int j = 0; j < crltns[i].size(); j++) {
 609             for (int o = 0; o < crltns[i][j].size(); o++) {
 610                 if (std::fabs(crltns[i][j][o] - crltns_test[i][j][o]) > 0) {
 611                     std::cout << "error " << i << " " << j << " " << o << " " << std::fabs(crltns[i][j][o] - crltns_test[i][j][o]) << "\n";
 612                 } else {
 613                     ola001++;
 614                 }
 615             }
 616         }
 617     }
 618     std::cout << "ola001 = " << ola001 << "\n";
 619     */
 620
 621     graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, reference, crltns, static_cast< double >(crl_border), static_cast< double >(eff_rad), m, k);
 622     std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
 623
 624     graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
 625
 626     std::cout << "routs evaluation: end\n";
 627
 628     make_rout_file(static_cast< double >(crl_border), index, rout_new, (OutPutName + "_routs.txt").c_str());
 629     std::cout << "corelation routs printed\n";
 630
 631     make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
 632     std::cout << "corelation routs' pairs' graphics printed\n";
 633
 634     /*std::cout << "extra params\n";
 635     std::vector< double > diffusion;
 636     evaluate_diffusion(trajectory, diffusion);
 637     make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);*/
 638
 639     std::cout << "Finish Analysis - end\n\n";
 640 }
 641
 642 void
 643 SpaceTimeCorr::writeOutput()
 644 {
 645
 646 }
 647
 648 /*! \brief
 649  * The main function for the analysis template.
 650  */
 651 int
 652 main(int argc, char *argv[])
 653 {
 654     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
 655 }