src/fitng.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 <chrono>
  45 #include <string.h>
  46 #include <algorithm>
  47 #include <cstdio>
  48
  49 #include <omp.h>
  50
  51 #include <gromacs/trajectoryanalysis.h>
  52 #include <gromacs/pbcutil/pbc.h>
  53 #include <gromacs/utility/smalloc.h>
  54 #include <gromacs/math/do_fit.h>
  55 #include <gromacs/fileio/trxio.h>
  56 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  57
  58 using namespace gmx;
  59
  60 using gmx::RVec;
  61
  62 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) {
  63     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) +
  64                     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) +
  65                     pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
  66 }
  67
  68 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) {
  69     double t01, t02, t03, t04, t05, t06, t07;
  70     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);
  71     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);
  72     t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
  73     t04 = sqrt(t01 + t02 + t03);
  74     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);
  75     t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
  76     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);
  77
  78     //#pragma omp parallel sections
  79     //{
  80         //#pragma omp section
  81             F += t04;
  82         //#pragma omp section
  83             Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
  84         //#pragma omp section
  85             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);
  86         //#pragma omp section
  87             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);
  88         //#pragma omp section
  89             Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
  90                 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 +
  91                 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);
  92         //#pragma omp section
  93             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 +
  94                 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 -
  95                 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
  96         //#pragma omp section
  97             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 -
  98                 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);
  99     //}
 100 }
 101
 102 /*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) {
 103     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;
 104     AmLfa = A - L * fa;
 105     BmLfb = B - L * fb;
 106     CmLfc = C - L * fc;
 107     BIXpXmLfx = bix + x - L * fx;
 108     BIYpYmLfy = biy + y - L * fy;
 109     BIZpZmLfz = biz + z - L * fz;
 110     t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
 111     t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
 112     t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
 113     t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
 114     t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
 115     t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
 116     t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
 117
 118     t8 = fc * sin(AmLfa) * sin(CmLfc);
 119     t9 = fa * cos(AmLfa) * cos(CmLfc);
 120     t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
 121     t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
 122     t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
 123
 124     t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
 125     t14 = fc * cos(AmLfa) * sin(CmLfc);
 126     t15 = fa * cos(CmLfc) * sin(AmLfa);
 127     t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
 128     t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
 129
 130     t18 = fx * cos(BmLfb) * sin(CmLfc);
 131     t19 = fc * cos(BmLfb) * cos(CmLfc);
 132     t20 = fb * sin(BmLfb) * sin(CmLfc);
 133
 134     fLd +=  (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 135                 (   BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
 136                 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 137                 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
 138                     fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 139                 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 140                 (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)) +
 141                     BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
 142                     fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
 143                     fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
 144             (2 * sqrt(  t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
 145
 146
 147     fLdd += (   (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 148                 (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) +
 149                     2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
 150                     pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
 151                     2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
 152                 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
 153                     fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
 154                     BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
 155                     fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
 156                     cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
 157                     fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 158                 (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
 159                     fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
 160                     fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
 161                     t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 162                     2 * fy * t3 + 2 * fz * t4 +
 163                     2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 164                 (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
 165                     cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 166                 (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
 167                     pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
 168                     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) +
 169                     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) -
 170                     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) +
 171                     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) +
 172                     2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
 173                     fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
 174                     2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
 175                     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 +
 176                     pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 177                 (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
 178                 (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) +
 179                 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 180                     fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
 181                 (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 +
 182                     2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 183                 (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 184                 (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
 185                     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) +
 186                     pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
 187                     2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
 188                     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) -
 189                     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) -
 190                     2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
 191                     pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
 192             (2 * sqrt(  t5 + t6 + t7)) -
 193             (   (2 *    (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 194                         (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
 195                         2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 196                         (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 197                             fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
 198                         2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 199                         (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)) +
 200                         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)) -
 201                         fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
 202                 (   (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
 203                     (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)) +
 204                         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)) -
 205                         fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
 206                     (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
 207                     (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) +
 208                         t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
 209                     (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
 210                         fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
 211             (4 * pow    (   t5 + t6 + t7, 3 / 2));
 212
 213 }*/
 214
 215 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
 216     #pragma omp parallel
 217     {
 218         #pragma omp for schedule(dynamic)
 219         for (int i = 0; i < b.size(); i++) {
 220             RVec temp = b[i];
 221             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);
 222             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);
 223             b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
 224         }
 225     }
 226 }
 227
 228 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
 229     midA[0] = 0;
 230     midA[1] = 0;
 231     midA[2] = 0;
 232
 233     midB[0] = 0;
 234     midB[1] = 0;
 235     midB[2] = 0;
 236
 237     for (int i = 0; i < pairs.size(); i++) {
 238         rvec_inc(midA, a[pairs[i].first]);
 239         rvec_inc(midB, b[pairs[i].second]);
 240     }
 241     midA[0] /= pairs.size();
 242     midA[1] /= pairs.size();
 243     midA[2] /= pairs.size();
 244
 245     midB[0] /= pairs.size();
 246     midB[1] /= pairs.size();
 247     midB[2] /= pairs.size();
 248 }
 249
 250 /*float*/double FrameMidLength(std::vector< RVec > a) {
 251     RVec b;
 252     b[0] = 0;
 253     b[1] = 0;
 254     b[2] = 0;
 255     for (int i = 0; i < a.size(); i++) {
 256         b += a[i];
 257     }
 258     return b.norm();
 259 }
 260
 261 void reset_coord(std::vector< RVec > &frame, RVec move) {
 262     for (int i = 0; i < frame.size(); i++) {
 263         frame[i] -= move;
 264     }
 265 }
 266
 267 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs, double FtCnst) {
 268     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
 269     RVec ma, mb;
 270     CalcMid(a, b, ma, mb, pairs);
 271
 272     reset_coord(a, ma);
 273     reset_coord(b, mb);
 274
 275     rvec_dec(ma, mb);
 276
 277     double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
 278     for (int i = 0; i < pairs.size(); i++) {
 279         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);
 280     }
 281     if (FtCnst == 0) {
 282         FtCnst = 0.0001;
 283     }
 284     if (f1 == 0) {
 285         return;
 286     } else {
 287         //int count = 0;
 288         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
 289             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 290             for (int i = 0; i < pairs.size(); i++) {
 291                 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);
 292             }
 293             while (true) {
 294                 f2 = 0;
 295                 for (int i = 0; i < pairs.size(); i++) {
 296                     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);
 297                 }
 298                 //count++;
 299                 if (f2 < f1) {
 300                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 301                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 302                     break;
 303                 } else {
 304                     l *= 0.99; // may be put somewhere l > 0
 305                 }
 306             }
 307         }
 308         //std::cout << " " << count << "\n";
 309         ApplyFit(b, x, y, z, A, B, C);
 310     }
 311 }
 312
 313 void printPDBtraj(const char* Fname, std::vector< std::vector < RVec > > trj, std::vector< std::vector< std::pair< int, int > > > prs, std::vector< int > ndx) {
 314     FILE* a = fopen (Fname, "a");
 315     for (int i = 1; i < trj.size(); i++) {
 316         fprintf(a, "TITLE     test t= %9.5f step= %5d\nMODEL %d\n", (/*float*/double)(i + 1), i + 1, i + 1);
 317         for (int j = 0; j < prs.size(); j++) {
 318             for (int k = 0; k < prs[j].size(); k++) {
 319                 fprintf(a, "ATOM  %5d  CA   CA %5d    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s  \n", ndx[prs[j][k].first], ndx[prs[j][k].first], (trj[i][prs[j][k].first][0] * 10), (trj[i][prs[j][k].first][1] * 10), (trj[i][prs[j][k].first][2] * 10), 1.0, 20.0, "    ", " ");
 320             }
 321         }
 322         fprintf(a, "TER\nENDMDL\n");
 323     }
 324     //fprintf(a, "END\n");
 325     fclose(a);
 326 }
 327
 328 bool myfunction (const std::vector< std::pair< int, int > > i, const std::vector< std::pair< int, int > > j) {
 329     return i.size() < j.size();
 330 }
 331
 332 void domain_reading(SelectionList sList, std::vector< int > ind, std::vector< std::vector < std::pair< int, int > > > &d_ind_num) {
 333     d_ind_num.resize(sList.size() - 1);
 334     for (int i = 1; i < sList.size(); i++) {
 335         d_ind_num[i - 1].resize(0);
 336         for (ArrayRef< const int >::iterator ai = sList[i].atomIndices().begin(); (ai < sList[i].atomIndices().end()); ai++) {
 337             d_ind_num[i - 1].push_back(std::make_pair(*ai, 0));
 338         }
 339     }
 340
 341     for (int i = 0; i < d_ind_num.size(); i++) {
 342         int k;
 343         for (int j = 0; j < d_ind_num[i].size(); j++) {
 344             k = 0;
 345             while (ind[k] != d_ind_num[i][j].first) {
 346                 k++;
 347             }
 348             d_ind_num[i][j].second = k;
 349         }
 350     }
 351
 352     // its a point of interest as we can either take biggest domains of the given set and work from those or we take them in given order
 353     //std::sort(d_ind_num.begin(), d_ind_num.end(), myfunction);
 354
 355     for (int i = 0; i < d_ind_num.size(); i++) {
 356         if (d_ind_num[i].size() >= 4) {
 357             for (int k = 0; k < d_ind_num[i].size(); k++) {
 358                 for (int j = i + 1; j < d_ind_num.size(); j++) {
 359                     for (int x = 0; x < d_ind_num[j].size(); x++) {
 360                         if (d_ind_num[j][x].first == d_ind_num[i][k].first) {
 361                             d_ind_num[j].erase(d_ind_num[j].begin() + x);
 362                             x--;
 363                         }
 364                     }
 365                 }
 366             }
 367         }
 368     }
 369 }
 370
 371 void domain_expansion(std::vector< int > indx, std::vector< std::vector < std::pair< int, int > > > &d_ind_num, std::vector< RVec > frame) {
 372     std::vector< bool > flag;
 373     flag.resize(indx.size(), 1);
 374     for (int i = 0; i < d_ind_num.size(); i++) {
 375         for (int j = 0; j < d_ind_num[i].size(); j++) {
 376             flag[d_ind_num[i][j].second] = 0;
 377         }
 378     }
 379     /*float*/double dist;
 380     int x3 = -1;
 381     for (int i = 0; i < indx.size(); i++) {
 382         if (flag[i]) {
 383             dist = 0;
 384             for (int x1 = 0; x1 < d_ind_num.size(); x1++) {
 385                 for (int x2 = 0; x2 < d_ind_num[x1].size(); x2++) {
 386                     if (dist == 0) {
 387                         dist = (frame[i] - frame[d_ind_num[0][0].first]).norm();
 388                         x3 = 0;
 389                     }
 390                     if ((frame[i] - frame[d_ind_num[x1][x2].first]).norm() < dist) {
 391                         x3 = x1;
 392                         dist = (frame[i] - frame[d_ind_num[x1][x2].first]).norm();
 393                     }
 394                 }
 395             }
 396             d_ind_num[x3].push_back(std::make_pair(indx[i], i));
 397         }
 398     }
 399 }
 400
 401 void MakeDomainFitPairs(std::vector< std::vector< std::pair< int, int > > > &p, std::vector< std::vector< std::pair< int, int > > > d) {
 402     p.resize(d.size());
 403     for (int i = 0; i < d.size(); i++) {
 404         p[i].resize(0);
 405         for (int j = 0; j < d[i].size(); j++) {
 406             p[i].push_back(std::make_pair(d[i][j].second, d[i][j].second));
 407         }
 408     }
 409 }
 410
 411 void TrajFitting(std::vector< std::vector < RVec > > &trj, double FC, std::vector< std::vector< std::pair< int, int > > > prs) {
 412
 413     // its needed to be thought through on whats shared and whats private
 414     //#pragma omp parallel for shared(trj, prs, FC)/* private(clone)*/
 415     for (int i = 1; i < trj.size()/*2*/; i++) {
 416         std::vector< std::vector < RVec > > clone;
 417         clone.resize(0);
 418         clone.resize(prs.size(), trj[i]);
 419         //std::chrono::time_point<std::chrono::system_clock> now1, now2;
 420         for (int j = 0; j < prs.size(); j++) {
 421             //now1 = std::chrono::system_clock::now();
 422             MyFitNew(trj[0], clone[j], prs[j], FC);
 423             //now2 = std::chrono::system_clock::now();
 424             //std::cout << " " << std::chrono::duration_cast<std::chrono::milliseconds>(now2 - now1).count() << "\n";
 425             for (int k = 0; k < prs[j].size(); k++) {
 426                 trj[i][prs[j][k].first] = clone[j][prs[j][k].first];
 427             }
 428         }
 429         clone.resize(0);
 430         //std::cout << "frame " << i << " fitted\n";
 431     }
 432 }
 433
 434 class Fitng : public TrajectoryAnalysisModule
 435 {
 436     public:
 437
 438         Fitng();
 439         virtual ~Fitng();
 440
 441         //! Set the options and setting
 442         virtual void initOptions(IOptionsContainer          *options,
 443                                  TrajectoryAnalysisSettings *settings);
 444
 445         //! First routine called by the analysis framework
 446         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 447         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 448                                   const TopologyInformation        &top);
 449
 450         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 451                                          const t_trxframe                   &fr);
 452
 453         //! Call for each frame of the trajectory
 454         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 455         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 456                                   TrajectoryAnalysisModuleData *pdata);
 457
 458         //! Last routine called by the analysis framework
 459         // virtual void finishAnalysis(t_pbc *pbc);
 460         virtual void finishAnalysis(int nframes);
 461
 462         //! Routine to write output, that is additional over the built-in
 463         virtual void writeOutput();
 464
 465     private:
 466         SelectionList                                               sel_;
 467
 468         std::vector< int >                                          index;
 469         const int                                                   *ind;
 470         std::vector< std::vector< std::pair< int, int > > >         domains_index_number;
 471         std::vector< std::vector < RVec > >                         trajectory;
 472         std::vector< RVec >                                         reference;
 473
 474         std::vector< t_trxframe >                                   saved_traj;
 475         //Selection                                                   selec;
 476         std::vector< std::vector< std::pair< int, int > > >         pairs;
 477         int                                                         frames = 0;
 478         int                                                         method = 1;
 479         double                                                      FitConst = 0;
 480         std::string                                                 OutPutTrjName;
 481         t_trxstatus                                                 *op;
 482 };
 483
 484 Fitng::Fitng(): TrajectoryAnalysisModule()
 485 {
 486 }
 487
 488 Fitng::~Fitng()
 489 {
 490 }
 491
 492 void
 493 Fitng::initOptions(IOptionsContainer          *options,
 494                      TrajectoryAnalysisSettings *settings)
 495 {
 496     static const char *const desc[] = {
 497         "[THISMODULE] to be done"
 498     };
 499     // Add the descriptive text (program help text) to the options
 500     settings->setHelpText(desc);
 501     // Add option for selecting a subset of atoms
 502     //options->addOption(SelectionOption("select")
 503     //                       .store(&selec).required()
 504     //                       .description("Atoms that are considered as part of the excluded volume"));
 505     // Add option for selection list
 506     options->addOption(SelectionOption("select_residue_and_domains").storeVector(&sel_)
 507                            .required().dynamicMask().multiValue()
 508                            .description("Selected group and domains based on it"));
 509     // Add option for Fit constant
 510     options->addOption(DoubleOption("FitConst")
 511                             .store(&FitConst)
 512                             .description("Fitting untill diff <= FitConst, by default == 0.00001"));
 513     // Add option for output pdb traj (meh edition)
 514     options->addOption(StringOption("trj_out")
 515                             .store(&OutPutTrjName)
 516                             .description("transformed trajectory"));
 517     // Add option for trajectory Fit transformation
 518     options->addOption(IntegerOption("method")
 519                             .store(&method)
 520                             .description("1 - fitting all -> all | 2 - keep only domains, all -> all | 3 - keep only domains, domain -> domain | 4 - extend domains to full coverage, domain -> domain"));
 521     // Control input settings
 522     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 523     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 524     settings->setPBC(true);
 525 }
 526
 527 // /home/toluk/Develop/samples/reca_test_ground/
 528 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL0' -n '/home/toluk/Develop/samples/reca_test_ground/t0.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/FullCAFit.pdb' -e 1000
 529 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsFit.pdb' -e 1000 -method 2
 530 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsOnlyFit.pdb' -e 1000 -method 3
 531 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsWExtraOnlyFit.pdb' -e 1000 -method 4
 532
 533 // -s '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.xtc' -sf '/home/titov_ai/Develop/samples/reca_rd/SLCa' -n '/home/titov_ai/Develop/samples/reca_rd/CA.ndx' -pdb_out '/home/titov_ai/Develop/samples/reca_rd/test.xtc' -method 1
 534
 535 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/test.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/testfit.xtc' -method 1
 536
 537 // -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_trjconvfit_time.xtc' -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_prefit_time.xtc' -method 1 -e 500 -FitConst 0.000001
 538 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1
 539
 540 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
 541
 542
 543 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test1k.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
 544
 545
 546 void
 547 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 548 {
 549     index.resize(0);
 550     for (ArrayRef< const int >::iterator ai = sel_[0].atomIndices().begin(); (ai < sel_[0].atomIndices().end()); ai++) {
 551         index.push_back(*ai);
 552     }
 553
 554     reference.resize(0);
 555     if (top.hasFullTopology()) {
 556         for (int i = 0; i < index.size(); i++) {
 557             reference.push_back(top.x().at(index[i]));
 558         }
 559     }
 560     //std::cout << reference.size() << "\n important \n";
 561
 562     saved_traj.resize(0);
 563
 564     if (sel_.size() > 1 && method >= 2) {
 565         domain_reading(sel_, index, domains_index_number);
 566     }
 567 }
 568
 569 void
 570 Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
 571 {
 572     trajectory.resize(0);
 573     std::vector < RVec > a;
 574     a.resize(0);
 575     trajectory.push_back(reference);
 576     /*for (int i = 0; i < index.size(); i++) {
 577         trajectory[0].push_back(fr.x[index[i]]);
 578     }*/
 579     trajectory.push_back(a);
 580     switch (method) {
 581         case 1:
 582             pairs.resize(1);
 583             for (int i = 0; i < index.size(); i++) {
 584                 pairs[0].push_back(std::make_pair(i, i));
 585             }
 586             break;
 587         case 2:
 588             pairs.resize(1);
 589             for (int i = 0; i < domains_index_number.size(); i++) {
 590                 for (int j = 0; j < domains_index_number[i].size(); j++) {
 591                     pairs[0].push_back(std::make_pair(domains_index_number[i][j].second, domains_index_number[i][j].second));
 592                 }
 593             }
 594             break;
 595         case 3:
 596             MakeDomainFitPairs(pairs, domains_index_number);
 597             break;
 598         case 4:
 599             domain_expansion(index, domains_index_number, trajectory[0]);
 600             MakeDomainFitPairs(pairs, domains_index_number);
 601             break;
 602     }
 603
 604     op = open_trx(OutPutTrjName.c_str(), "w+");
 605 }
 606
 607 void
 608 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 609                       TrajectoryAnalysisModuleData *pdata)
 610 {
 611     trajectory[1].resize(0);
 612     for (int i = 0; i < index.size(); i++) {
 613         trajectory[1].push_back(fr.x[index[i]]);
 614     }
 615     TrajFitting(trajectory, FitConst, pairs);
 616
 617     t_trxframe c = fr;
 618     for (int i = 0; i < index.size(); i++) {
 619         copy_rvec(trajectory[1][i].as_vec(), c.x[index[i]]);
 620     }
 621     //trajectory.push_back(trajectory[1]);
 622     if (write_trxframe_indexed(op, static_cast<const t_trxframe *>(&c), index.size(), static_cast<const int *>(&index.front())/*saved_traj[i].index*/, NULL) != 0) {
 623         std::cout << "\nFileOutPut error\nframe number = " << frnr << "\n";
 624     }
 625
 626     //printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
 627 }
 628
 629
 630 void
 631 Fitng::finishAnalysis(int nframes)
 632 {
 633     std::cout << "\nGame Over\n";
 634 }
 635
 636 void
 637 Fitng::writeOutput()
 638 {
 639     //printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
 640 }
 641
 642
 643 /*! \brief
 644  * The main function for the analysis template.
 645  */
 646 int
 647 main(int argc, char *argv[])
 648 {
 649     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
 650 }