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