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