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 #include <cfloat>
  53
  54 #include <gromacs/trajectoryanalysis.h>
  55 #include <gromacs/utility/smalloc.h>
  56 #include <gromacs/math/do_fit.h>
  57 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  58 #include "gromacs/topology/topology.h"
  59 #include "gromacs/topology/index.h"
  60
  61 using namespace gmx;
  62
  63 using gmx::RVec;
  64
  65 double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
  66     return  sqrt(   pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
  67                     pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
  68                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  69 }
  70
  71 void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
  72     double t01, t02, t03, t04, t05, t06, t07;
  73     t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
  74     t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
  75     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  76     t04 = sqrt(t01 + t02 + t03);
  77     t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
  78     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  79     t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
  80     F += t04;
  81     Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  82     Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
  83     Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
  84     Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  85             2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
  86             2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
  87     Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
  88             2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
  89             2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  90     Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
  91             2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
  92 }
  93
  94 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
  95     double t0 = 0, t1 = 0, t2 = 0;
  96     for (unsigned int i = 0; i < b.size(); i++) {
  97         t0 = static_cast< double >(b[i][0]);
  98         t1 = static_cast< double >(b[i][1]);
  99         t2 = static_cast< double >(b[i][2]);
 100         b[i][0] = static_cast< float >((t2 + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (t1 + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (t0 + x));
 101         b[i][1] = static_cast< float >((t1 + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (t2 + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (t0 + x));
 102         b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
 103     }
 104 }
 105
 106 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
 107     midA[0] = 0;
 108     midA[1] = 0;
 109     midA[2] = 0;
 110
 111     midB[0] = 0;
 112     midB[1] = 0;
 113     midB[2] = 0;
 114
 115     for (unsigned int i = 0; i < pairs.size(); i++) {
 116         //midA += a[pairs[i].first];
 117         //midB += b[pairs[i].second];
 118         rvec_inc(midA, a[pairs[i].first]);
 119         rvec_inc(midB, b[pairs[i].second]);
 120     }
 121     midA[0] /= pairs.size();
 122     midA[1] /= pairs.size();
 123     midA[2] /= pairs.size();
 124
 125     midB[0] /= pairs.size();
 126     midB[1] /= pairs.size();
 127     midB[2] /= pairs.size();
 128 }
 129
 130 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
 131     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 132     RVec ma, mb;
 133     CalcMid(a, b, ma, mb, pairs);
 134     //ma -= mb;
 135     rvec_dec(ma, mb);
 136     double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
 137     for (unsigned int i = 0; i < pairs.size(); i++) {
 138         f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 139                 static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 140     }
 141     if (FtCnst == 0) {
 142         FtCnst = 0.00001;
 143     }
 144     if (f1 == 0) {
 145         return;
 146     } else {
 147         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
 148             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
 149             for (unsigned int i = 0; i < pairs.size(); i++) {
 150                 searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
 151                                static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 152                                static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
 153                                static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
 154             }
 155             while (f2 != f1) {
 156                 f2 = 0;
 157                 for (unsigned int i = 0; i < pairs.size(); i++) {
 158                     f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
 159                             static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
 160                             static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
 161                 }
 162                 if (f2 < f1) {
 163                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 164                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 165                     break;
 166                 } else {
 167                     l *= 0.50;
 168                     if (l == DBL_MIN) { //DBL_TRUE_MIN
 169                         break;
 170                     }
 171                 }
 172             }
 173         }
 174         ApplyFit(b, x, y, z, A, B, C);
 175     }
 176 }
 177
 178 //const int tau_jump = 10;
 179
 180 const std::vector< std::vector< std::vector< double > > > read_correlation_matrix_file(const char* file_name, unsigned long size, unsigned int length)
 181 {
 182     FILE *file;
 183     std::vector< std::vector< std::vector< double > > > crrlts;
 184     file = std::fopen(file_name, "r+");
 185     std::vector< double > c;
 186     c.resize(size, 0);
 187     std::vector< std::vector< double > > b;
 188     b.resize(size, c);
 189     crrlts.resize(length + 1, b);
 190     //char a[100];
 191     std::cout << "reading correlations from a file:\n";
 192     for (unsigned int i = 0; i < length + 1; i++) {
 193         int t0 = 0, t1 = std::fscanf(file, "%d\n", &t0);
 194         std::cout << t0 << " | ";
 195         if (i % 100 == 0) {
 196             std::cout << "\n";
 197         }
 198         for (unsigned int j = 0; j < size; j++) {
 199             for (unsigned int f = 0; f < size; f++) {
 200                 int t2 = std::fscanf(file, "%lf ", &crrlts[i][j][f]);
 201             }
 202         }
 203     }
 204     std::fclose(file);
 205     std::cout << "\n";
 206     return crrlts;
 207 }
 208
 209 template< typename T >
 210 unsigned long posSrch(std::vector< T > a, T b) {
 211     for (unsigned i01 = 0; i01 < a.size(); i01++) {
 212         if (a[i01] == b) {
 213             return i01;
 214             break;
 215         }
 216     }
 217     return a.size() + 1;
 218 }
 219
 220 void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std::string > resNames,
 221                     std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout,
 222                     const char* file_name01, const unsigned int tauD, const unsigned int window_size, int t_jump)
 223 {
 224     FILE *file01;
 225     file01 = std::fopen(file_name01, "w+");
 226
 227     unsigned int f, pre = 0;
 228     std::vector< std::tuple< int, int, std::vector< int > > > table;
 229
 230     for (unsigned int i01 = 0; i01 < rout.size(); i01++) {
 231         f = i01;
 232         for (unsigned int i02 = i01 + 1; i02 <rout.size(); i02++) {
 233             if (rout[i02] == rout[i01]) {
 234                 f = i02;
 235             } else {
 236                 break;
 237             }
 238         }
 239         if (i01 == f) {
 240             std::fprintf(file01, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", i01 * window_size / t_jump, crl_border, tauD, window_size);
 241         } else {
 242             std::fprintf(file01, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", i01 * window_size / t_jump, f * window_size / t_jump, crl_border, tauD, window_size);
 243         }
 244         for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
 245             for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
 246                 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);
 247             }
 248             std::fprintf(file01, "\n");
 249         }
 250         table.resize(0);
 251         for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
 252             for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
 253                 bool flag1 = true, flag2 = true;
 254                 for (unsigned int i04 = 0; i04 < table.size(); i04++) {
 255                     if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].first]) {
 256                         std::get<1>(table[i04])++;
 257                         std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].second]);
 258                         flag1 = false;
 259                     }
 260                     if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].second]) {
 261                         std::get<1>(table[i04])++;
 262                         std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].first]);
 263                         flag2 = false;
 264                     }
 265                 }
 266                 if (flag1) {
 267                     std::vector< int > a;
 268                     a.resize(0);
 269                     a.push_back(indx[rout[i01][i02][i03].second]);
 270                     table.push_back(std::make_tuple(indx[rout[i01][i02][i03].first], 1, a));
 271                 }
 272                 if (flag2) {
 273                     std::vector< int > a;
 274                     a.resize(0);
 275                     a.push_back(indx[rout[i01][i02][i03].first]);
 276                     table.push_back(std::make_tuple(indx[rout[i01][i02][i03].second], 1, a));
 277                 }
 278             }
 279         }
 280
 281         for (unsigned int i02 = 0; i02 < table.size(); i02++) {
 282             std::fprintf(file01, "residue %s connections %d | ", (resNames[posSrch(indx, std::get<0>(table[i02]))]).c_str(), std::get<1>(table[i02]));
 283             for (unsigned int i03 = 0; i03 < std::get<2>(table[i02]).size(); i03++) {
 284                 std::fprintf(file01, "%s ", (resNames[posSrch(indx, std::get<2>(table[i02])[i03])]).c_str());
 285             }
 286             std::fprintf(file01, "\n");
 287         }
 288         if (pre != 0) {
 289             std::vector< std::vector< std::pair< unsigned int, unsigned int > > > temp01, temp02;
 290             std::vector< std::pair< unsigned int, unsigned int > > temp03, temp04, temp05;
 291             temp01 = rout[pre];
 292             temp02 = rout[f];
 293             temp03.resize(0);
 294             temp04.resize(0);
 295             for (unsigned int i02 = 0; i02 < temp01.size(); i02++) {
 296                 for (unsigned int i03 = 0; i03 < temp01[i02].size(); i03++) {
 297                     temp03.push_back(temp01[i02][i03]);
 298                 }
 299             }
 300             for (unsigned int i02 = 0; i02 < temp02.size(); i02++) {
 301                 for (unsigned int i03 = 0; i03 < temp02[i02].size(); i03++) {
 302                     temp04.push_back(temp02[i02][i03]);
 303                 }
 304             }
 305             std::sort(temp03.begin(), temp03.end());
 306             std::sort(temp04.begin(), temp04.end());
 307             temp05.resize(0);
 308             std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
 309             std::fprintf(file01, "minus:\n");
 310             for (unsigned int i02 = 0; i02 < temp05.size(); i02++) {
 311                 std::fprintf(file01, "cgo_arrow (id %3d), (id %3d), radius=0.075\n", indx[temp05[i02].first] + 1, indx[temp05[i02].second] + 1);
 312             }
 313             temp05.resize(0);
 314             std::set_difference(temp04.begin(), temp04.end(), temp03.begin(), temp03.end(), std::inserter(temp05, temp05.begin()));
 315             std::fprintf(file01, "plus:\n");
 316             for (unsigned int i02 = 0; i02 < temp05.size(); i02++) {
 317                 std::fprintf(file01, "cgo_arrow (id %3d), (id %3d), radius=0.075\n", indx[temp05[i02].first] + 1, indx[temp05[i02].second] + 1);
 318             }
 319         }
 320         pre = f;
 321         i01 = f;
 322     }
 323     std::fclose(file01);
 324 }
 325
 326 bool isitsubset (std::vector< int > a, std::vector< int > b) {
 327     if (b.size() == 0) {
 328         return true;
 329     }
 330     std::sort(a.begin(), a.end());
 331     std::sort(b.begin(), b.end());
 332     unsigned int k = 0;
 333     for (unsigned int i = 0; i < a.size(); i++) {
 334         if (b[k] == a[i]) {
 335             k++;
 336         }
 337     }
 338     if (k == b.size()) {
 339         return true;
 340     } else {
 341         return false;
 342     }
 343 }
 344
 345 float minDistToSubset(const std::vector< RVec > set, std::vector < unsigned int > subSet, int pos) {
 346     float a = 9000;
 347     for (unsigned int j = 0; j < subSet.size(); j++) {
 348         if (a > gmx::norm((set[pos]-set[j]).as_vec())) {
 349             a = gmx::norm((set[pos]-set[j]).as_vec());
 350         }
 351     }
 352     return a;
 353 }
 354
 355 void correlation_evaluation(const std::vector< RVec > ref,
 356                             const std::vector< std::vector< RVec > > trajectory,
 357                             const unsigned int tauE, const unsigned int window_size, const char* file_name, int t_jump,
 358                             std::vector< std::vector < std::vector < unsigned int > > > sls) {
 359     std::vector< std::vector< std::vector< double > > > crl;
 360     std::vector< std::vector< RVec > > trj;
 361
 362     for (unsigned int istart = 0; istart < trajectory.size() - window_size - tauE; istart += window_size / t_jump) {
 363         trj.clear();
 364         trj.resize(window_size + tauE);
 365         for (unsigned int i = 0; i < window_size + tauE; i++) {
 366             for (unsigned int j = 0; j < trajectory[i].size(); j++) {
 367                 trj[i].push_back(trajectory[i + istart][j]);
 368             }
 369         }
 370
 371         /*
 372          *
 373          * fitting block
 374          * we add spare atoms to existing domains based on proximity
 375          *
 376          */
 377
 378         std::vector< std::vector< std::pair< unsigned int, unsigned int > > > prs;
 379         int prsCount = 0;
 380         float prsTemp = 9000;
 381         for (unsigned int i0 = 0; i0 < ref.size(); i0++) {
 382             for (unsigned int i = 0; i < sls[istart / (window_size / t_jump)].size(); i++) {
 383                 if (prsTemp > minDistToSubset(ref, sls[istart / (window_size / t_jump)][i], i0)) {
 384                     prsTemp = minDistToSubset(ref, sls[istart / (window_size / t_jump)][i], i0);
 385                     prsCount = i;
 386                 }
 387                 if (prsTemp < 0.1) { // to avoid 0. error
 388                     prsCount = -1;
 389                     break;
 390                 }
 391             }
 392             if (prsCount != -1) {
 393                 sls[istart / (window_size / t_jump)][prsCount].push_back(i0);
 394             }
 395         }
 396
 397         prs.resize(sls[istart / (window_size / t_jump)].size());
 398         for (unsigned int i = 0; i < sls[istart / (window_size / t_jump)].size(); i++) {
 399             for (unsigned int j = 0; j < sls[istart / (window_size / t_jump)][i].size(); j++) {
 400                 prs[i].push_back(std::make_pair(sls[istart / (window_size / t_jump)][i][j], sls[istart / (window_size / t_jump)][i][j]));
 401             }
 402         }
 403
 404         std::vector< std::vector< std::vector< RVec > > > trjTemp;
 405         trjTemp.resize(prs.size());
 406         for (unsigned int i = 0; i < prs.size(); i++) {
 407             trjTemp[i] = trj;
 408             for (unsigned int j = 0; j < trjTemp[i].size(); j++) {
 409                 MyFitNew(ref, trjTemp[i][j], prs[i], 0);
 410                 for (unsigned int k = 0; k < prs[i].size(); k++) {
 411                     trj[j][prs[i][k].first] = trjTemp[i][j][prs[i][k].first];
 412                 }
 413             }
 414         }
 415
 416         /*
 417          *
 418          * fitting block
 419          *
 420          */
 421
 422
 423         crl.resize(tauE + 1);
 424         for (unsigned int i = 0; i < crl.size(); i++) {
 425             crl[i].resize(trajectory.front().size());
 426             for (unsigned int j = 0; j < trajectory.front().size(); j++) {
 427                 crl[i][j].resize(trajectory.front().size());
 428                 for (unsigned int k = 0; k < trajectory.front().size(); k++) {
 429                     crl[i][j][k] = 0.;
 430                 }
 431             }
 432         }
 433         #pragma omp parallel for ordered schedule(dynamic) shared(crl) firstprivate(trj, ref)
 434         for (unsigned int i = 0; i <= tauE; i++) {
 435             std::vector< std::vector< double > > a, b, c;
 436             std::vector< double > d;
 437             d.resize(0);
 438             a.resize(0);
 439             b.resize(0);
 440             c.resize(0);
 441             for (unsigned int j = 0; j < trajectory.front().size(); j++) {
 442                 a.push_back(d);
 443                 b.push_back(d);
 444                 c.push_back(d);
 445                 for (unsigned int k = 0; k < trajectory.front().size(); k++) {
 446                     a[j].push_back(0.);
 447                     b[j].push_back(0.);
 448                     c[j].push_back(0.);
 449                 }
 450             }
 451             for (unsigned int j = 0; j < trj.size() - tauE; j++) {
 452                 for (unsigned int f1 = 0; f1 < trj.front().size(); f1++) {
 453                     for (unsigned int f2 = 0; f2 < trj.front().size(); f2++) {
 454                         RVec temp1, temp2;
 455                         temp1 = trj[j][f1] - ref[f1];
 456                         temp2 = trj[j + i][f2] - ref[f2];
 457                         a[f1][f2] +=    (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
 458                                         static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
 459                                         static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
 460                         b[f1][f2] +=    (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
 461                                         static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
 462                                         static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
 463                         c[f1][f2] +=    (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
 464                                         static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
 465                                         static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
 466                     }
 467                 }
 468             }
 469             for (unsigned int j = 0; j < trj.front().size(); j++) {
 470                 for (unsigned int f = 0; f < trj.front().size(); f++) {
 471                     crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 472                 }
 473             }
 474         }
 475         #pragma omp barrier
 476
 477         for (unsigned int i1 = 0; i1 < crl.size(); i1++) {
 478             for (unsigned int i2 = 0; i2 < crl[i1].size(); i2++) {
 479                 for (unsigned int i3 = 0; i3 < crl[i1][i2].size(); i3++) {
 480                     crl[i1][i2][i3] = std::round(crl[i1][i2][i3] * 10000) / 10000;
 481                 }
 482             }
 483         }
 484
 485         FILE *file;
 486         file = std::fopen(file_name, "a");
 487         for (unsigned int i = 0; i < crl.size(); i++) {
 488             std::fprintf(file, "%d %d\n", istart, i);
 489             for (unsigned int j = 0; j < crl[i].size(); j++) {
 490                 for (unsigned int f = 0; f < crl[i][j].size(); f++) {
 491                     std::fprintf(file, "%.6f ", crl[i][j][f]); //~16
 492                 }
 493                 std::fprintf(file, "\n");
 494             }
 495         }
 496         std::fclose(file);
 497         if (istart % 1000 == 0 && istart > 0) {
 498             std::cout << "\n";
 499         }
 500         std::cout << " | done " << istart << " | ";
 501     }
 502     std::cout << "\n";
 503 }
 504
 505 void graph_calculation( std::vector< std::vector< std::pair< double, long int > > > &graph, std::vector< std::vector< unsigned int > > &s_graph,
 506                         std::vector< std::vector< std::pair< unsigned int, unsigned int > > > &s_graph_rbr,
 507                         const std::vector< RVec > ref,
 508                         const std::vector< std::vector< std::vector<double > > > crl,
 509                         const double crl_b, const double e_rad, const unsigned int tauB, const unsigned int tauE) {
 510     graph.resize(ref.size());
 511     for (unsigned int i = 0; i < ref.size(); i++) {
 512         graph[i].resize(ref.size(), std::make_pair(0, -1));
 513     }
 514     RVec temp;
 515     for (unsigned int i = tauB; i <= tauE; i++) {
 516         for (unsigned int j = 0; j < ref.size(); j++) {
 517             for (unsigned int f = j; f < ref.size(); f++) {
 518                 temp = ref[i] - ref[f];
 519                 if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b &&
 520                     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]))) {
 521                     if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
 522                         graph[j][f].first = crl[i][j][f];
 523                     } else {
 524                         graph[j][f].first = crl[i][f][j];
 525                     }
 526                     graph[j][f].second = i;
 527                 }
 528             }
 529         }
 530     }
 531     std::vector< bool > graph_flags;
 532     graph_flags.resize(ref.size(), true);
 533     std::vector< unsigned int > a;
 534     a.resize(0);
 535     std::vector< std::pair< unsigned int, unsigned int > > b;
 536     b.resize(0);
 537     std::vector< unsigned int > que1, que2, que3;
 538     for (unsigned int i = 0; i < ref.size(); i++) {
 539         if (graph_flags[i]) {
 540             s_graph.push_back(a);
 541             s_graph_rbr.push_back(b);
 542             que1.resize(0);
 543             que2.resize(0);
 544             que3.resize(0);
 545             que1.push_back(i);
 546             que3.push_back(i);
 547             graph_flags[i] = false;
 548             while(que1.size() > 0) {
 549                 que2.clear();
 550                 for (unsigned int k = 0; k < que1.size(); k++) {
 551                     for (unsigned int j = 0; j < ref.size(); j++) {
 552                         if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
 553                             que2.push_back(j);
 554                             graph_flags[j] = false;
 555                         }
 556                     }
 557                 }
 558                 que1 = que2;
 559                 for (unsigned int j = 0; j < que2.size(); j++) {
 560                     que3.push_back(que2[j]);
 561                 }
 562             }
 563             s_graph.back() = que3;
 564             for (unsigned int j = 0; j < que3.size(); j++) {
 565                 for (unsigned int k = 0; k < ref.size(); k++) {
 566                     if (graph[que3[j]][k].second > -1) {
 567                         s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
 568                     }
 569                 }
 570             }
 571         }
 572     }
 573 }
 574
 575 bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
 576     return i.second < j.second;
 577 }
 578
 579 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< unsigned int, unsigned int > > > &rout_n, const unsigned long indxSize,
 580                                 const std::vector< std::vector< std::pair< double, long > > > graph, const std::vector< std::vector< unsigned int > > s_graph,
 581                                 const std::vector< std::vector< std::pair< unsigned int, unsigned int > > > s_graph_rbr) {
 582     std::vector< double > key;
 583     std::vector< long > path;
 584     std::vector< std::pair< unsigned int, double > > que;
 585     std::vector< std::pair< unsigned int, unsigned int > > a;
 586     unsigned int v;
 587     for (unsigned int i = 0; i < s_graph.size(); i++) {
 588         key.resize(0);
 589         path.resize(0);
 590         que.resize(0);
 591         v = 0;
 592         if (s_graph[i].size() > 2) {
 593             key.resize(indxSize, 2);
 594             path.resize(indxSize, -1);
 595             key[s_graph[i][0]] = 0;
 596             for (unsigned int j = 0; j < s_graph[i].size(); j++) {
 597                 que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
 598             }
 599             std::sort(que.begin(), que.end(), myfunction);
 600             while (!que.empty()) {
 601                 v = que[0].first;
 602                 que.erase(que.begin());
 603                 for (unsigned int j = 0; j < s_graph_rbr[i].size(); j++) {
 604                     long u = -1;
 605                     if (s_graph_rbr[i][j].first == v) {
 606                         u = s_graph_rbr[i][j].second;
 607                     } else if (s_graph_rbr[i][j].second == v) {
 608                         u = s_graph_rbr[i][j].first;
 609                     }
 610                     bool flag = false;
 611                     unsigned long pos = 0;
 612                     for (unsigned long k = 0; k < que.size(); k++) {
 613                         if (que[k].first == u) {
 614                             flag = true;
 615                             pos = k;
 616                             k = que.size() + 1;
 617                         }
 618                     }
 619                     if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
 620                         path[static_cast< unsigned long >(u)] = v;
 621                         key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
 622                         que[pos].second = key[static_cast< unsigned long >(u)];
 623                         sort(que.begin(), que.end(), myfunction);
 624                     }
 625                 }
 626             }
 627             a.resize(0);
 628             rout_n.push_back(a);
 629             for (unsigned int j = 0; j < indxSize; j++) {
 630                 if (path[j] != -1) {
 631                     rout_n.back().push_back(std::make_pair(j, path[j]));
 632                 }
 633             }
 634         }
 635     }
 636 }
 637
 638 gmx::RVec evaluate_com(std::vector< RVec > frame) {
 639     RVec temp;
 640     temp[0] = 0;
 641     temp[1] = 0;
 642     temp[2] = 0;
 643     for (unsigned int i = 0; i < frame.size(); i++) {
 644         temp += frame[i];
 645     }
 646     temp[0] /= frame.size();
 647     temp[1] /= frame.size();
 648     temp[2] /= frame.size();
 649     return temp;
 650 }
 651
 652 /*! \brief
 653  * \ingroup module_trajectoryanalysis
 654  */
 655
 656 class SpaceTimeCorr : public TrajectoryAnalysisModule
 657 {
 658     public:
 659
 660         SpaceTimeCorr();
 661         virtual ~SpaceTimeCorr();
 662
 663         //! Set the options and setting
 664         virtual void initOptions(IOptionsContainer          *options,
 665                                  TrajectoryAnalysisSettings *settings);
 666
 667         //! First routine called by the analysis framework
 668         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 669         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 670                                   const TopologyInformation        &top);
 671
 672         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 673                                          const t_trxframe                   &fr);
 674
 675         //! Call for each frame of the trajectory
 676         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 677         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 678                                   TrajectoryAnalysisModuleData *pdata);
 679
 680         //! Last routine called by the analysis framework
 681         // virtual void finishAnalysis(t_pbc *pbc);
 682         virtual void finishAnalysis(int nframes);
 683
 684         //! Routine to write output, that is additional over the built-in
 685         virtual void writeOutput();
 686
 687     private:
 688
 689         SelectionList                                               sel_;
 690
 691         std::vector< std::vector< RVec > >                          trajectory;
 692         std::vector< RVec >                                         reference;
 693
 694         std::vector< std::vector < std::vector < unsigned int > > > sels;
 695
 696         std::vector< int >                                          index;
 697         std::vector< std::string >                                  resNms;
 698
 699         int                                                         frames              = 0;
 700         int                                                         tau                 = 0;    // selectable
 701
 702         int                                                         tau_jump            = 10;
 703
 704         int                                                         window              = 0;
 705         float                                                       crl_border          = 0;    // selectable
 706         float                                                       eff_rad             = 1.5;  // selectable
 707         std::string                                                 OutPutName;                 // selectable
 708         int                                                         mode                = 0;    // selectable
 709         // Copy and assign disallowed by base.
 710 };
 711
 712 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
 713 {
 714 }
 715
 716 SpaceTimeCorr::~SpaceTimeCorr()
 717 {
 718 }
 719
 720 void
 721 SpaceTimeCorr::initOptions(IOptionsContainer          *options,
 722                      TrajectoryAnalysisSettings *settings)
 723 {
 724     static const char *const desc[] = {
 725         "[THISMODULE] to be done"
 726     };
 727     // Add the descriptive text (program help text) to the options
 728     settings->setHelpText(desc);
 729     // Add option for output file name
 730     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 731     //                        .store(&fnNdx_).defaultBasename("domains")
 732     //                        .description("Index file from the domains"));
 733     // Add option for working mode
 734     options->addOption(gmx::IntegerOption("mode")
 735                             .store(&mode)
 736                             .description("default 0 - creation of sigle matrixes' file | 1 - reading prepared correlation file to make graphs and etc."));
 737     // Add option for tau constant
 738     options->addOption(gmx::IntegerOption("tau")
 739                             .store(&tau)
 740                             .description("max number of frames to shift to evaluate correlations | all starts from array [0 ; tau]"));
 741     // Add option for window constant
 742     options->addOption(gmx::IntegerOption("window")
 743                             .store(&window)
 744                             .description("size of window in frames that correlation evaluation based on"));
 745     // Add option for crl_border constant
 746     options->addOption(FloatOption("crl")
 747                             .store(&crl_border)
 748                             .description("make graphs based on corrs > crl const"));
 749     // Add option for eff_rad constant
 750     options->addOption(FloatOption("ef_rad")
 751                             .store(&eff_rad)
 752                             .description("effective radius for atoms to see each other to evaluate corrs"));
 753     // Add option for selection list
 754     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 755                            .required().dynamicMask().multiValue()
 756                            .description("Domains to form rigid skeleton"));
 757     // Add option for output file names
 758     options->addOption(StringOption("out_put")
 759                             .store(&OutPutName)
 760                             .description("<your name here> + <local file tag>.txt"));
 761     // Control input settings
 762     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 763     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 764     settings->setPBC(true);
 765 }
 766
 767 void
 768 SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 769 {
 770     ArrayRef< const int > atomind = sel_.front().atomIndices();
 771     index.resize(0);
 772     for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 773         index.push_back(*ai);
 774         resNms.push_back(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name) + std::to_string((top.atoms()->atom[*ai].resind + 1)));
 775     }
 776     trajectory.resize(0);
 777
 778     std::string str(sel_.back().name());
 779     unsigned int tempSize = std::stoul(str.substr(13, 6)) / (window / tau_jump);
 780
 781     sels.resize(tempSize);
 782
 783     for (unsigned int i = 1; i < sel_.size(); i++) {
 784         std::string str(sel_[i].name());
 785         atomind = sel_[i].atomIndices();
 786         std::vector< unsigned int > a;
 787         a.resize(0);
 788         sels[std::stoul(str.substr(13, 6)) / (window / tau_jump)].push_back(a);
 789         for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 790             sels[std::stoul(str.substr(13, 6)) / (window / tau_jump)].back().push_back(*ai);
 791         }
 792     }
 793
 794     reference.resize(0);
 795     if (top.hasFullTopology()) {
 796         for (unsigned int i = 0; i < index.size(); i++) {
 797             reference.push_back(top.x().at(index[i]));
 798         }
 799     }
 800 }
 801
 802 void
 803 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
 804 {
 805 }
 806
 807 // -s *.tpr -f *.xtc -n *.ndx -sf 'name' -mode <> -tau <> -window <> -crl <> -ef_rad <> -out_put <>
 808
 809 void
 810 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 811 {
 812     trajectory.resize(trajectory.size() + 1);
 813     trajectory.back().resize(index.size());
 814     for (unsigned int i = 0; i < index.size(); i++) {
 815         trajectory.back()[i] = fr.x[index[i]];
 816     }
 817     frames++;
 818 }
 819
 820 void
 821 SpaceTimeCorr::finishAnalysis(int nframes)
 822 {
 823     std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout_new;
 824     unsigned int k = 500, m = 0;
 825     if (tau > 0) {
 826         k = static_cast< unsigned int>(tau);
 827     }
 828     if (window == 0) {
 829         window = static_cast< int >(k * 2);
 830     }
 831
 832     // need to double check the formula
 833     //std::cout << "\nProgram requires " << 4*index.size()*index.size()*static_cast< unsigned int>(tau)*(trajectory.size() - static_cast< unsigned int>(window) - k)/tau_jump/1024/1024/1024 << " Gb of free space (may be x2)";
 834
 835     if (mode == 0) {
 836         std::cout << "\nCorrelation's evaluation - start\n";
 837
 838         correlation_evaluation(reference, trajectory, k, static_cast< unsigned int >(window), (OutPutName + "-DumpData.txt").c_str(), tau_jump, sels);
 839
 840         std::cout << "Correlation's evaluation - end\n";
 841     } else if (mode == 1) {
 842         rout_new.resize((trajectory.size() - static_cast< unsigned long >(window) - k)/tau_jump);
 843         FILE *file;
 844         file = std::fopen((OutPutName + "-DumpData.txt").c_str(), "r+");
 845         //#pragma omp parallel for ordered schedule(dynamic) shared(rout_new) firstprivate(reference, crl_border, eff_rad, m, k)
 846         for(unsigned long i = 0; i < trajectory.size() - static_cast< unsigned int>(window) - k; i+= static_cast< unsigned int >(window) / tau_jump) {
 847             std::vector< std::vector< std::vector< double > > > crltns;
 848             std::vector< std::vector< std::pair< double, long int > > > graph;
 849             std::vector< std::vector< unsigned int > > sub_graph;
 850             std::vector< std::vector< std::pair< unsigned int, unsigned int > > > sub_graph_rbr;
 851
 852             //#pragma omp ordered
 853             //{
 854                 crltns.resize(k + 1);
 855                 for (unsigned int i1 = 0; i1 < crltns.size(); i1++) {
 856                     crltns[i1].resize(trajectory.front().size());
 857                     for (unsigned int i2 = 0; i2 < trajectory.front().size(); i2++) {
 858                         crltns[i1][i2].resize(trajectory.front().size());
 859                         for (unsigned int i3 = 0; i3 < trajectory.front().size(); i3++) {
 860                             crltns[i1][i2][i3] = 0.;
 861                         }
 862                     }
 863                 }
 864                 for (unsigned int i1 = 0; i1 < k + 1; i1++) {
 865                     int t0, t1, t2 = std::fscanf(file, "%d %d\n", &t0, &t1);
 866                     for (unsigned int i2 = 0; i2 < trajectory.front().size(); i2++) {
 867                         for (unsigned int i3 = 0; i3 < trajectory.front().size(); i3++) {
 868                             t2 = std::fscanf(file, "%lf ", &crltns[i1][i2][i3]);
 869                         }
 870                     }
 871                 }
 872                 if (i % 10 == 0 && i > 0) {
 873                     std::cout << "\n";
 874                 }
 875                 std::cout << " | mtrxs " << i << " red | ";
 876             //}
 877
 878             graph_calculation(graph, sub_graph, sub_graph_rbr, reference, crltns, static_cast< double >(crl_border), static_cast< double >(eff_rad), m, k);
 879             graph_back_bone_evaluation(rout_new[i / (static_cast< unsigned int >(window) / tau_jump)], reference.size(), graph, sub_graph, sub_graph_rbr);
 880
 881         }
 882         //#pragma omp barrier
 883         std::cout << "\n Almost - end\n";
 884         make_rout_file(static_cast< double >(crl_border), index, resNms, rout_new, (OutPutName + "-routs.txt").c_str(), static_cast< unsigned int >(tau), static_cast< unsigned int >(window), tau_jump);
 885     }
 886     std::cout << "Finish Analysis - end\n\n";
 887 }
 888
 889 void
 890 SpaceTimeCorr::writeOutput()
 891 {
 892
 893 }
 894
 895 /*! \brief
 896  * The main function for the analysis template.
 897  */
 898 int
 899 main(int argc, char *argv[])
 900 {
 901     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
 902 }