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