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