src/spacetimecorr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43 #include <iostream>
  44 #include <algorithm>
  45 #include <omp.h>
  46 #include <set>
  47
  48 #include <gromacs/trajectoryanalysis.h>
  49 #include <gromacs/utility/smalloc.h>
  50 #include <gromacs/math/do_fit.h>
  51
  52 #include <vector>
  53 #include <math.h>
  54
  55 //#include "fittingn.cpp"
  56
  57 using namespace gmx;
  58
  59 using gmx::RVec;
  60
  61 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) {
  62     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) +
  63                     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) +
  64                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  65 }
  66
  67 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) {
  68     double t01, t02, t03, t04, t05, t06, t07;
  69     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);
  70     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);
  71     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  72     t04 = sqrt(t01 + t02 + t03);
  73     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);
  74     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  75     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);
  76
  77     //#pragma omp parallel sections
  78     //{
  79         //#pragma omp section
  80             F += t04;
  81         //#pragma omp section
  82             Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  83         //#pragma omp section
  84             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);
  85         //#pragma omp section
  86             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);
  87         //#pragma omp section
  88             Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  89                 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 +
  90                 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);
  91         //#pragma omp section
  92             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 +
  93                 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 -
  94                 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  95         //#pragma omp section
  96             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 -
  97                 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);
  98     //}
  99 }
 100
 101 void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
 102     double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
 103     AmLfa = A - L * fa;
 104     BmLfb = B - L * fb;
 105     CmLfc = C - L * fc;
 106     BIXpXmLfx = bix + x - L * fx;
 107     BIYpYmLfy = biy + y - L * fy;
 108     BIZpZmLfz = biz + z - L * fz;
 109     t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
 110     t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
 111     t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
 112     t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
 113     t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
 114     t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
 115     t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
 116
 117     t8 = fc * sin(AmLfa) * sin(CmLfc);
 118     t9 = fa * cos(AmLfa) * cos(CmLfc);
 119     t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
 120     t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
 121     t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
 122
 123     t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
 124     t14 = fc * cos(AmLfa) * sin(CmLfc);
 125     t15 = fa * cos(CmLfc) * sin(AmLfa);
 126     t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
 127     t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
 128
 129     t18 = fx * cos(BmLfb) * sin(CmLfc);
 130     t19 = fc * cos(BmLfb) * cos(CmLfc);
 131     t20 = fb * sin(BmLfb) * sin(CmLfc);
 132
 133     fLd +=  (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 134                 (   BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
 135                 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 136                 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
 137                     fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 138                 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 139                 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
 140                     BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
 141                     fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
 142                     fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
 143             (2 * sqrt(  t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
 144
 145
 146     fLdd += (   (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 147                 (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
 148                     2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
 149                     pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
 150                     2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
 151                 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
 152                     fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
 153                     BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
 154                     fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
 155                     cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
 156                     fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 157                 (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
 158                     fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
 159                     fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
 160                     t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 161                     2 * fy * t3 + 2 * fz * t4 +
 162                     2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 163                 (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
 164                     cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 165                 (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
 166                     pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
 167                     pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
 168                     2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
 169                     pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
 170                     pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
 171                     2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
 172                     fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 173                     2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
 174                     fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
 175                     pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 176                 (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
 177                 (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
 178                 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 179                     fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
 180                 (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 181                     2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 182                 (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 183                 (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
 184                     pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
 185                     pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
 186                     2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
 187                     pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
 188                     2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
 189                     2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
 190                     pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
 191             (2 * sqrt(  t5 + t6 + t7)) -
 192             (   (2 *    (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 193                         (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
 194                         2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 195                         (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 196                             fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 197                         2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 198                         (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
 199                         BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 200                         fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
 201                 (   (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 202                     (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
 203                         BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 204                         fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 205                     (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 206                     (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
 207                         t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 208                     (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 209                         fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
 210             (4 * pow    (   t5 + t6 + t7, 3 / 2));
 211
 212 }
 213
 214
 215 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
 216     RVec temp;
 217     //#pragma omp parallel
 218     //{
 219         //#pragma omp for schedule(dynamic)
 220         for (int i = 0; i < b.size(); i++) {
 221             temp = b[i];
 222             b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
 223             b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
 224             b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
 225         }
 226     //}
 227 }
 228
 229 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
 230     midA[0] = 0;
 231     midA[1] = 0;
 232     midA[2] = 0;
 233
 234     midB[0] = 0;
 235     midB[1] = 0;
 236     midB[2] = 0;
 237
 238     for (int i = 0; i < pairs.size(); i++) {
 239         rvec_inc(midA, a[pairs[i].first]);
 240         rvec_inc(midB, b[pairs[i].second]);
 241     }
 242     midA[0] /= pairs.size();
 243     midA[1] /= pairs.size();
 244     midA[2] /= pairs.size();
 245
 246     midB[0] /= pairs.size();
 247     midB[1] /= pairs.size();
 248     midB[2] /= pairs.size();
 249 }
 250
 251 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs) {
 252     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 253     RVec ma, mb;
 254     CalcMid(a, b, ma, mb, pairs);
 255     rvec_dec(ma, mb);
 256     double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
 257     for (int i = 0; i < pairs.size(); i++) {
 258         f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
 259     }
 260     if (f1 == 0) {
 261         return;
 262     } else {
 263         double fLd, fLdd;
 264         //l = 0;
 265         int count = 0;
 266             while (f1 - f2 < 0 || f1 - f2 > 0.00001) {
 267                 f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 268                 for (int i = 0; i < pairs.size(); i++) {
 269                     searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
 270                 }
 271                 /*do {
 272                     fLd = 0;
 273                     fLdd = 0;
 274                     for (int i = 0; i < pairs.size(); i++) {
 275                         searchL(fLd, fLdd, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0], b[pairs[i].second][1], b[pairs[i].second][2], x, y, z, A, B, C, l, fx, fy, fz, fa, fb, fc);
 276                     }
 277                     l -= fLd / fLdd;
 278                 } while (std::abs(fLd) > 1);*/
 279                 //std::cout << " Ltempo= " << Ltempo << " iter= " << Lco << "\n";
 280                 while (true) {
 281                     f2 = 0;
 282                     for (int i = 0; i < pairs.size(); i++) {
 283                         f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
 284                     }
 285                     count++;
 286                     if (f2 < f1) {
 287                         x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 288                         fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 289                         break;
 290                     } else {
 291                         l *= 0.85;
 292                     }
 293                 }
 294                 count++;
 295             }
 296         //std::cout << "iteration count: " << count << "\n";
 297         ApplyFit(b, x, y, z, A, B, C);
 298     }
 299 }
 300
 301 void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
 302 {
 303     FILE *file;
 304     file = std::fopen(file_name, "w+");
 305     for (int i = 1; i < trajectory.size(); i++) {
 306         std::fprintf(file, "MODEL%9d\n", i);
 307         for (int j = 0; j < trajectory[i].size(); j++) {
 308             std::fprintf(file, "ATOM  %5d   CA LYS     1    %8.3f%8.3f%8.3f  1.00  0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10);
 309         }
 310         std::fprintf(file, "ENDMDL\n");
 311     }
 312 }
 313
 314 void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
 315 {
 316     FILE *file;
 317     file = std::fopen(file_name, "w+");
 318     for (int i = start; i < correlations.size(); i++) {
 319         if (correlations.size() - start > 1) {
 320             std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
 321         }
 322         if (i % 50 == 0) {
 323             std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
 324         }
 325         for (int j = 0; j < correlations[i].size(); j++) {
 326             for (int f = 0; f < correlations[i][j].size(); f++) {
 327                 std::fprintf(file, "%3.2f ", correlations[i][j][f]);
 328             }
 329             std::fprintf(file, "\n");
 330         }
 331     }
 332     std::fclose(file);
 333 }
 334
 335 void make_rout_file2(float crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)
 336 {
 337     FILE *file;
 338     file = std::fopen(file_name, "w+");
 339     std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
 340     for (int i = 0; i < rout.size(); i++) {
 341         for (int j = 0; j < rout[i].size(); j++) {
 342             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);
 343         }
 344         std::fprintf(file, "\n\n");
 345     }
 346     std::fclose(file);
 347 }
 348
 349 void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > > correlations,
 350                               std::vector< std::vector< std::pair< int, int > > > rout_pairs,
 351                               std::vector< int > indx,
 352                               const char* file_name)
 353 {
 354     FILE *file;
 355     file = std::fopen(file_name, "w+");
 356     for (int i = 0; i < rout_pairs.size(); i++) {
 357         for (int j = 0; j < rout_pairs[i].size(); j++) {
 358             std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first] + 1, indx[rout_pairs[i][j].second] + 1);
 359             for (int k = 0; k < correlations.size(); k++) {
 360                 std::fprintf(file, "%3.5f ", correlations[k][rout_pairs[i][j].first][rout_pairs[i][j].second]);
 361             }
 362             std::fprintf(file, "\n");
 363         }
 364     }
 365     std::fclose(file);
 366 }
 367
 368 bool mysortfunc (std::vector< int > a, std::vector< int > b) { return (a.size() > b.size()); }
 369
 370 bool isitsubset (std::vector< int > a, std::vector< int > b) {
 371     if (b.size() == 0) {
 372         return true;
 373     }
 374     std::sort(a.begin(), a.end());
 375     std::sort(b.begin(), b.end());
 376     int k = 0;
 377     for (int i = 0; i < a.size(); i++) {
 378         if (b[k] == a[i]) {
 379             k++;
 380         }
 381     }
 382     if (k == b.size()) {
 383         return true;
 384     } else {
 385         return false;
 386     }
 387 }
 388
 389 void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {
 390     crl.resize(tauE + 1);
 391     for (int i = 0; i < crl.size(); i++) {
 392         crl[i].resize(traj.front().size());
 393         for (int j = 0; j < crl[i].size(); j++) {
 394             crl[i][j].resize(traj.front().size(), 0);
 395         }
 396     }
 397     rvec temp_zero;
 398     temp_zero[0] = 0;
 399     temp_zero[1] = 0;
 400     temp_zero[2] = 0;
 401
 402     std::vector< float > d;
 403     d.resize(traj.front().size(), 0);
 404
 405     #pragma omp parallel shared(crl, temp_zero, d)
 406     {
 407         #pragma omp for schedule(dynamic)
 408         for (int i = tauS; i <= tauE; i += 1) {
 409             rvec *medx, *medy, temp1, temp2;
 410             real tmp;
 411             snew(medx, traj.front().size());
 412             snew(medy, traj.front().size());
 413             for (int i = 0; i < traj.front().size(); i++) {
 414                 copy_rvec(temp_zero, medx[i]);
 415                 copy_rvec(temp_zero, medy[i]);
 416             }
 417             for (int j = 0; j < traj.size() - i - 1; j++) {
 418                 for (int f = 0; f < traj.front().size(); f++) {
 419                     rvec_sub(traj[b_frame][f], traj[j][f], temp1);
 420                     rvec_inc(medx[f], temp1);
 421
 422                     rvec_sub(traj[b_frame][f], traj[j + i][f], temp1);
 423                     rvec_inc(medy[f], temp1);
 424                 }
 425             }
 426             for (int j = 0; j < traj.front().size(); j++) {
 427                 tmp = traj.size() - 1;
 428                 tmp -= i;
 429                 tmp = 1 / tmp;
 430                 //tmp = 1 / (traj.size() - 1 - i);
 431                 copy_rvec(medx[j], temp1);
 432                 svmul(tmp, temp1, medx[j]);
 433                 copy_rvec(medy[j], temp1);
 434                 svmul(tmp, temp1, medy[j]);
 435             }
 436             std::vector< std::vector< float > > a, b, c;
 437             a.resize(traj.front().size(), d);
 438             b.resize(traj.front().size(), d);
 439             c.resize(traj.front().size(), d);
 440             for (int j = 0; j < traj.size() - i - 1; j++) {
 441                 for (int f1 = 0; f1 < traj.front().size(); f1++) {
 442                     for (int f2 = 0; f2 < traj.front().size(); f2++) {
 443                         rvec_sub(traj[b_frame][f1], traj[j][f1], temp1);
 444                         rvec_dec(temp1, medx[f1]);
 445
 446                         rvec_sub(traj[b_frame][f2], traj[j + i][f2], temp2);
 447                         rvec_dec(temp2, medy[f2]);
 448
 449                         a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
 450                         b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
 451                         c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
 452                     }
 453                 }
 454             }
 455             for (int j = 0; j < traj.front().size(); j++) {
 456                 for (int f = 0; f < traj.front().size(); f++) {
 457                     crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 458                 }
 459             }
 460             sfree(medx);
 461             sfree(medy);
 462             std::cout << i << " corr done\n";
 463         }
 464     }
 465     #pragma omp barrier
 466 }
 467
 468 void domain_chopping(SelectionList seList, std::vector< int > indx, std::vector< std::vector< int > > &filled_domains, std::vector< RVec > frame) {
 469     ConstArrayRef< int > atomInd  = seList[0].atomIndices();
 470     std::vector< std::vector< int > > init_domains;
 471     init_domains.resize(seList.size());
 472     for (int i = 0; i < seList.size(); i++) {
 473         if (seList.size() != 1 && i == 0) {
 474             continue;
 475         }
 476         atomInd  = seList[i].atomIndices();
 477         for (ConstArrayRef<int>::iterator ai = atomInd.begin(); (ai < atomInd.end()); ai++) {
 478             init_domains[i].push_back(*ai);
 479         }
 480     }
 481     for (int i = 0; i < init_domains.size(); i++) {
 482         for (int j = 0; j < init_domains[i].size(); j++) {
 483             int k = 0;
 484             while (indx[k] != init_domains[i][j]) {
 485                 k++;
 486             }
 487             init_domains[i][j] = k;
 488         }
 489     }
 490     for (int i = 0; i < init_domains.size(); i++) {
 491         if (init_domains[i].size() >= 2) {
 492             filled_domains.push_back(init_domains[i]);
 493             for (int k = 0; k < init_domains[i].size(); k++) {
 494                 for (int j = i + 1; j < init_domains.size(); j++) {
 495                     for (int x = 0; x < init_domains[j].size(); x++) {
 496                         if (init_domains[j][x] == init_domains[i][k]) {
 497                             init_domains[j].erase(init_domains[j].begin() + x);
 498                         }
 499                     }
 500                 }
 501             }
 502         }
 503     }
 504     init_domains.resize(0);
 505     init_domains = filled_domains;
 506     std::vector< bool > flags;
 507     flags.resize(indx.size(), true);
 508     for (int i = 0; i < init_domains.size(); i++) {
 509         for (int j = 0; j < init_domains[i].size(); j++) {
 510             flags[init_domains[i][j]] = false;
 511         }
 512     }
 513     int a;
 514     rvec temp;
 515     for (int i = 0; i < indx.size(); i++) {
 516         if (flags[i]) {
 517             float dist = 90000001;
 518             for (int j = 0; j < init_domains.size(); j++) {
 519                 for (int k = 0; k < init_domains[j].size(); k++) {
 520                     rvec_sub(frame[i], frame[init_domains[j][k]], temp);
 521                     if (norm(temp) <= dist) {
 522                         dist = norm(temp);
 523                         a = j;
 524                     }
 525                 }
 526             }
 527             filled_domains[a].push_back(i);
 528             flags[i] = false;
 529         }
 530     }
 531 }
 532
 533 void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
 534                       std::vector< std::vector< RVec > > traj, int b_frame,
 535                       std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) {
 536     graph.resize(traj.front().size());
 537     for (int i = 0; i < traj.front().size(); i++) {
 538         graph[i].resize(traj.front().size(), std::make_pair(0, -1));
 539     }
 540     rvec temp;
 541     for (int i = 1; i <= tauE; i++) {
 542         for (int j = 0; j < traj.front().size(); j++) {
 543             for (int f = j; f < traj.front().size(); f++) {
 544                 copy_rvec(traj[b_frame][j], temp);
 545                 rvec_dec(temp, traj[b_frame][f]);
 546                 if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b && norm(temp) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
 547                     if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
 548                         graph[j][f].first = crl[i][j][f];
 549                     } else {
 550                         graph[j][f].first = crl[i][f][j];
 551                     }
 552                     graph[j][f].second = i;
 553                 }
 554             }
 555         }
 556     }
 557     std::cout << "crl analysed\n";
 558     std::vector< bool > graph_flags;
 559     graph_flags.resize(traj.front().size(), true);
 560     std::vector< int > a;
 561     a.resize(0);
 562     std::vector< std::pair< int, int > > b;
 563     b.resize(0);
 564     std::vector< int > que1, que2, que3;
 565     for (int i = 0; i < traj.front().size(); i++) {
 566         if (graph_flags[i]) {
 567             s_graph.push_back(a);
 568             s_graph_rbr.push_back(b);
 569             que1.resize(0);
 570             que2.resize(0);
 571             que3.resize(0);
 572             que1.push_back(i);
 573             que3.push_back(i);
 574             graph_flags[i] = false;
 575             while(que1.size() > 0) {
 576                 que2.clear();
 577                 for (int k = 0; k < que1.size(); k++) {
 578                     for (int j = 0; j < traj.front().size(); j++) {
 579                         if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
 580                             que2.push_back(j);
 581                             graph_flags[j] = false;
 582                         }
 583                     }
 584                 }
 585                 que1 = que2;
 586                 for (int j = 0; j < que2.size(); j++) {
 587                     que3.push_back(que2[j]);
 588                 }
 589             }
 590             s_graph.back() = que3;
 591             for (int j = 0; j < que3.size(); j++) {
 592                 for (int k = 0; k < traj.front().size(); k++) {
 593                     if (graph[que3[j]][k].second > -1) {
 594                         s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
 595                     }
 596                 }
 597             }
 598             //std::cout << s_graph.back().size() << "   ";
 599         }
 600     }
 601 }
 602
 603 bool myfunction (const std::pair< int, float > i, const std::pair< int, float > j) {
 604     return i.second < j.second;
 605 }
 606
 607 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,
 608                                 std::vector< std::vector< std::pair< float, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
 609     std::vector< float > key;
 610     std::vector< int > path;
 611     std::vector< std::pair< int, float > > que;
 612     std::vector< std::pair< int, int > > a;
 613     int v;
 614     for (int i = 0; i < s_graph.size(); i++) {
 615         key.resize(0);
 616         path.resize(0);
 617         que.resize(0);
 618         v = 0;
 619         if (s_graph[i].size() > 2) {
 620             key.resize(indxSize, 2);
 621             path.resize(indxSize, -1);
 622             key[s_graph[i][0]] = 0;
 623             for (int j = 0; j < s_graph[i].size(); j++) {
 624                 que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
 625             }
 626             std::sort(que.begin(), que.end(), myfunction);
 627             while (!que.empty()) {
 628                 v = que[0].first;
 629                 que.erase(que.begin());
 630                 for (int j = 0; j < s_graph_rbr[i].size(); j++) {
 631                     int u = -1;
 632                     if (s_graph_rbr[i][j].first == v) {
 633                         u = s_graph_rbr[i][j].second;
 634                     } else if (s_graph_rbr[i][j].second == v) {
 635                         u = s_graph_rbr[i][j].first;
 636                     }
 637                     bool flag = false;
 638                     int pos = -1;
 639                     for (int k = 0; k < que.size(); k++) {
 640                         if (que[k].first == u) {
 641                             flag = true;
 642                             pos = k;
 643                             k = que.size() + 1;
 644                         }
 645                     }
 646                     if (flag && key[u] > 1 - std::abs(graph[v][u].first)) {
 647                         path[u] = v;
 648                         key[u] = 1 - std::abs(graph[v][u].first);
 649                         que[pos].second = key[u];
 650                         sort(que.begin(), que.end(), myfunction);
 651                     }
 652                 }
 653             }
 654             a.resize(0);
 655             rout_n.push_back(a);
 656             for (int j = 0; j < indxSize; j++) {
 657                 if (path[j] != -1) {
 658                     rout_n.back().push_back(std::make_pair(j, path[j]));
 659                 }
 660             }
 661         }
 662     }
 663 }
 664
 665 /*! \brief
 666  * \ingroup module_trajectoryanalysis
 667  */
 668
 669 class Domains : public TrajectoryAnalysisModule
 670 {
 671     public:
 672
 673         Domains();
 674         virtual ~Domains();
 675
 676         //! Set the options and setting
 677         virtual void initOptions(IOptionsContainer          *options,
 678                                  TrajectoryAnalysisSettings *settings);
 679
 680         //! First routine called by the analysis framework
 681         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 682         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 683                                   const TopologyInformation        &top);
 684
 685         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 686                                          const t_trxframe                   &fr);
 687
 688         //! Call for each frame of the trajectory
 689         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 690         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 691                                   TrajectoryAnalysisModuleData *pdata);
 692
 693         //! Last routine called by the analysis framework
 694         // virtual void finishAnalysis(t_pbc *pbc);
 695         virtual void finishAnalysis(int nframes);
 696
 697         //! Routine to write output, that is additional over the built-in
 698         virtual void writeOutput();
 699
 700     private:
 701
 702         std::string                                                 fnNdx_;
 703         SelectionList                                               sel_;
 704
 705
 706         std::vector< std::vector< int > >                           domains_local;
 707         std::vector< std::vector< RVec > >                          trajectory;
 708         std::vector< std::vector< RVec > >                          frankenstein_trajectory;
 709
 710         std::vector< int >                                          index;
 711         std::vector< int >                                          numbers;
 712         int                                                         frames              = 0;
 713         int                                                         basic_frame         = 0;
 714         int                                                         tau                 = 0;
 715         float                                                      crl_border          = 0;
 716         float                                                      eff_rad             = 1.5;
 717         real                                                        **w_rls;
 718         std::vector< std::vector< std::pair< int, int > > >         w_rls2;
 719
 720         int                                                         domains_ePBC;
 721         // Copy and assign disallowed by base.
 722 };
 723
 724 Domains::Domains(): TrajectoryAnalysisModule()
 725 {
 726 }
 727
 728 Domains::~Domains()
 729 {
 730 }
 731
 732 void
 733 Domains::initOptions(IOptionsContainer          *options,
 734                      TrajectoryAnalysisSettings *settings)
 735 {
 736     static const char *const desc[] = {
 737         "[THISMODULE] to be done"
 738     };
 739     // Add the descriptive text (program help text) to the options
 740     settings->setHelpText(desc);
 741     // Add option for output file name
 742     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 743                             .store(&fnNdx_).defaultBasename("domains")
 744                             .description("Index file from the domains"));
 745     // Add option for tau constant
 746     options->addOption(gmx::IntegerOption("tau")
 747                             .store(&tau)
 748                             .description("number of frames for time travel"));
 749     // Add option for crl_border constant
 750     options->addOption(FloatOption("crl")
 751                             .store(&crl_border)
 752                             .description("make graph based on corrs > constant"));
 753     // Add option for eff_rad constant
 754     options->addOption(FloatOption("ef_rad")
 755                             .store(&eff_rad)
 756                             .description("effective radius for atoms to evaluate corrs"));
 757     // Add option for selection list
 758     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 759                            .required().dynamicMask().multiValue()
 760                            .description("Domains to form rigid skeleton"));
 761     // Control input settings
 762     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 763     settings->setPBC(true);
 764 }
 765
 766 void
 767 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 768                       const TopologyInformation        &top)
 769 {
 770     domains_ePBC = top.ePBC();
 771 }
 772
 773 void
 774 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 775                              const t_trxframe                       &fr)
 776 {
 777     //std::vector< std::vector< int > > domains;
 778     ConstArrayRef< int > atomind  = sel_[0].atomIndices();
 779     index.resize(0);
 780     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 781         index.push_back(*ai);
 782     }
 783     trajectory.resize(1);
 784     trajectory.back().resize(index.size());
 785
 786     for (int i = 0; i < index.size(); i++) {
 787         trajectory.back()[i] = fr.x[index[i]];
 788     }
 789
 790     domain_chopping(sel_, index, domains_local, trajectory.back());
 791
 792     frankenstein_trajectory.resize(trajectory.size());
 793     frankenstein_trajectory.back() = trajectory.back();
 794     trajectory.resize(trajectory.size() + 1);
 795     trajectory.back().resize(0);
 796     trajectory.back().resize(index.size());
 797
 798     w_rls2.resize(domains_local.size() + 1);
 799     for (int i = 0; i < domains_local.size(); i++) {
 800         w_rls2[i].resize(domains_local[i].size());
 801         for (int j = 0; j < domains_local[i].size(); j++) {
 802             w_rls2[i][j] = std::make_pair(domains_local[i][j], domains_local[i][j]);
 803         }
 804     }
 805     w_rls2[domains_local.size()].resize(index.size());
 806     for (int i = 0; i < index.size(); i++) {
 807         w_rls2.back()[i] = std::make_pair(i, i);
 808     }
 809     frames++;
 810 }
 811
 812 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5'  -tau 5 -crl 0.10
 813 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/CorrsTestDomainsNfit.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionListDomainsNFit' -tau 5000 -crl 0.75 -ef_rad 1
 814 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/TestCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 100 -crl 0.75 -ef_rad 1
 815
 816 void
 817 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 818                       TrajectoryAnalysisModuleData *pdata)
 819 {
 820     for (int i = 0; i < index.size(); i++) {
 821         trajectory.back()[i] = fr.x[index[i]];
 822     }
 823     frankenstein_trajectory.resize(frankenstein_trajectory.size() + 1);
 824     frankenstein_trajectory.back().resize(index.size());
 825     #pragma omp parallel shared(frankenstein_trajectory, trajectory)
 826     {
 827         #pragma omp for schedule(dynamic)
 828         for (int i = 0; i < domains_local.size(); i++) {
 829             std::vector< RVec > traj = trajectory.back();
 830             MyFitNew(trajectory[basic_frame], traj, w_rls2[i]);
 831             for (int j = 0; j < domains_local[i].size(); j++) {
 832                 frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
 833             }
 834         }
 835     }
 836     #pragma omp barrier
 837     std::cout << "frame № " << frames + 1 <<" analysed\n";
 838     frames++;
 839 }
 840
 841 void
 842 Domains::finishAnalysis(int nframes)
 843 {
 844     std::vector< std::vector< std::vector< float > > > crltns;
 845     int k = 1000, m = 0;
 846     if (tau > 0) {
 847         k = tau;
 848     }
 849
 850     /*
 851      *
 852      *
 853      *
 854     */
 855
 856     std::cout << "Correlation's evaluation - start\n\n";
 857     correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);
 858     make_correlation_file(crltns, "full crltns file Ca.txt", 0);
 859     std::cout << "Correlation's evaluation - end\n\n";
 860
 861     /*
 862      *
 863      *
 864      *
 865     */
 866
 867     //number of corelations in the matrixes
 868     /*for (int i1 = 0; i1 < 100; i1++) {
 869         int i5 = 0;
 870         for (int i2 = 1; i2 < crltns.size(); i2++) {
 871             for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
 872                 for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {
 873                     if (std::abs(crltns[i2][i3][i4]) >= (float)i1 / 100 && i3 != i4) {
 874                         i5++;
 875                     }
 876                 }
 877             }
 878         }
 879         std::cout << i1 << " - " << i5 << " | ";
 880         if (i1 % 10 == 0) {
 881             std::cout << "\n" ;
 882         }
 883     }*/
 884
 885     std::cout << "graph evaluation: start\n";
 886
 887     std::vector< std::vector< std::pair< float, int > > > graph;
 888     std::vector< std::vector< int > > sub_graph;
 889     std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr;
 890
 891
 892     graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, k);
 893
 894
 895     std::vector< std::vector < std::pair< int, int > > > rout_new;
 896
 897     std::cout << "graph evaluation: end\n";
 898     std::cout << "routs evaluation: start\n";
 899
 900     graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
 901
 902     std::cout << "routs evaluation: end\n";
 903
 904
 905     // Routs_DomainsNFit_5000_0.75_1_t.txt
 906     make_rout_file2(crl_border, index, rout_new, "Routs_Ca.txt");
 907     make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Ca.txt");
 908
 909     std::cout << "Finish Analysis - end\n\n";
 910 }
 911
 912 void
 913 Domains::writeOutput()
 914 {
 915
 916 }
 917
 918 /*! \brief
 919  * The main function for the analysis template.
 920  */
 921 int
 922 main(int argc, char *argv[])
 923 {
 924     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 925 }