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