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 <omp.h>
  46 #include <thread>
  47 #include <string>
  48 #include <math.h>
  49 #include <algorithm>
  50
  51 #include <gromacs/trajectoryanalysis.h>
  52 #include <gromacs/pbcutil/pbc.h>
  53 #include <gromacs/pbcutil/rmpbc.h>
  54 #include <gromacs/utility/smalloc.h>
  55 #include <gromacs/math/vectypes.h>
  56 #include <gromacs/math/vec.h>
  57 #include <gromacs/math/do_fit.h>
  58
  59 using namespace gmx;
  60
  61 using gmx::RVec;
  62
  63 double F (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
  64     return  sqrt(   pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
  65                     pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2) +
  66                     pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2)  );
  67 }
  68
  69 double Fx (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
  70     return  -(2 * cos(B) * cos(C) * (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
  71             2 * sin(B) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y)) +
  72             2 * cos(B) * sin(C) * (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x))) /
  73             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
  74             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
  75             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
  76 }
  77
  78 double Fy (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
  79     return  -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) *
  80             (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) -
  81             2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) *
  82             (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) +
  83             2 * cos(B) * sin(A) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
  84             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
  85             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
  86             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
  87 }
  88
  89 double Fz (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
  90     return  -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) *
  91             (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
  92             2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) *
  93             (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
  94             2 * cos(A) * cos(B) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
  95             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
  96             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
  97             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
  98 }
  99
 100 double FA (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
 101     return  -(2 * (cos(A) * cos(B) * (biy + y) - cos(B) * sin(A) * (biz + z)) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y)) -
 102             2 * ((biy + y) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + (biz + z) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) *
 103             (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
 104             2 * ((biy + y) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + (biz + z) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) *
 105             (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x))) /
 106             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
 107             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
 108             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
 109 }
 110
 111 double FB (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
 112     return  -(2 * (cos(A) * cos(B) * sin(C) * (biz + z) - sin(B) * sin(C) * (bix + x) + cos(B) * sin(A) * sin(C) * (biy + y)) *
 113             (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
 114             2 * (cos(A) * cos(B) * cos(C) * (biz + z) - cos(C) * sin(B) * (bix + x) + cos(B) * cos(C) * sin(A) * (biy + y)) *
 115             (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
 116             2 * (cos(B) * (bix + x) + sin(A) * sin(B) * (biy + y) + cos(A) * sin(B) * (biz + z)) *
 117             (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
 118             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
 119             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
 120             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
 121 }
 122
 123 double FC (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
 124     return  (2 * ((biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (bix + x)) *
 125             (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
 126             2 * ((biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (bix + x)) *
 127             (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x))) /
 128             (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
 129             pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
 130             pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
 131 }
 132
 133 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
 134     RVec temp;
 135     for (int i = 0; i < b.size(); i++) {
 136         temp = b[i];
 137         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);
 138         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);
 139         b[i][2] = cos(A)*cos(B)*(temp[2] + z) - sin(B)*(temp[0] + x) + cos(B)*sin(A)*(temp[1] + y);
 140     }
 141 }
 142
 143 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
 144     midA[0] = 0;
 145     midA[1] = 0;
 146     midA[2] = 0;
 147
 148     midB[0] = 0;
 149     midB[1] = 0;
 150     midB[2] = 0;
 151
 152     for (int i = 0; i < pairs.size(); i++) {
 153         rvec_inc(midA, a[pairs[i].first]);
 154         rvec_inc(midB, b[pairs[i].second]);
 155     }
 156     midA[0] /= pairs.size();
 157     midA[1] /= pairs.size();
 158     midA[2] /= pairs.size();
 159
 160     midB[0] /= pairs.size();
 161     midB[1] /= pairs.size();
 162     midB[2] /= pairs.size();
 163 }
 164
 165 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs) {
 166     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0;
 167     double f21 = 0, f22 = 0, f23 = 0, f24 = 0, f25 = 0, f26 = 0;
 168     double l1 = 1, l2 = 1, l3 = 1, l4 = 1, l5 = 1, l6 = 1, l = 1;
 169     RVec ma, mb;
 170     CalcMid(a, b, ma, mb, pairs);
 171     rvec_dec(ma, mb);
 172     double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
 173
 174     int count;
 175     while (true) {
 176         for (int i = 0; i < pairs.size(); i++) {
 177             f1 += F(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);
 178             fx += Fx(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);
 179             fy += Fy(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);
 180             fz += Fz(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);
 181             fa += FA(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);
 182             fb += FB(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);
 183             fc += FC(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);
 184         }
 185         while (true) {
 186             f2 = 0;
 187             for (int i = 0; i < pairs.size(); i++) {
 188                 f2 += F(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 - l * fx, y - l * fy, z - l * fz, A - l * fa, B - l * fb, C - l * fc);
 189                 /*
 190                 f2 += F(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 - l1 * fx, y - l2 * fy, z - l3 * fz, A - l4 * fa, B - l5 * fb, C - l6 * fc);
 191                 f21 += F(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 - l1 * fx, y, z, A, B, C);
 192                 f22 += F(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 - l2 * fy, z, A, B, C);
 193                 f23 += F(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 - l3 * fz, A, B, C);
 194                 f24 += F(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 - l4 * fa, B, C);
 195                 f25 += F(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 - l5 * fb, C);
 196                 f26 += F(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 - l6 * fc);
 197                 */
 198             }
 199             /*count++;
 200             if (count % 100 == 0) {
 201                 std::cout << f1 << " " << f2 << "\n";
 202             }*/
 203             if (f2 < f1) {
 204                 x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 205                 //x -= l1 * fx; y -= l2 * fy; z -= l3 * fz; A -= l4 * fa; B -= l5 * fb; C -= l6 * fc;
 206                 fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 207                 break;
 208             } else {
 209                 l *= 0.85;
 210                 /*
 211                 if (f21 > f1) {
 212                     l1 /= 2;
 213                 }
 214                 if (f22 > f1) {
 215                     l2 /= 2;
 216                 }
 217                 if (f23 > f1) {
 218                     l3 /= 2;
 219                 }
 220                 if (f24 > f1) {
 221                     l4 /= 2;
 222                 }
 223                 if (f25 > f1) {
 224                     l5 /= 2;
 225                 }
 226                 if (f26 > f1) {
 227                     l6 /= 2;
 228                 }
 229                 */
 230             }
 231         }
 232         if (f1 - f2 > 0 && f1 - f2 < 0.00001) {
 233             break;
 234         }
 235         f1 = 0; f2 = 0;
 236     }
 237     ApplyFit(b, x, y, z, A, B, C);
 238 }
 239
 240 class Fitng : public TrajectoryAnalysisModule
 241 {
 242     public:
 243
 244         Fitng();
 245         virtual ~Fitng();
 246
 247         //! Set the options and setting
 248         virtual void initOptions(IOptionsContainer          *options,
 249                                  TrajectoryAnalysisSettings *settings);
 250
 251         //! First routine called by the analysis framework
 252         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 253         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 254                                   const TopologyInformation        &top);
 255
 256         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 257                                          const t_trxframe                   &fr);
 258
 259         //! Call for each frame of the trajectory
 260         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 261         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 262                                   TrajectoryAnalysisModuleData *pdata);
 263
 264         //! Last routine called by the analysis framework
 265         // virtual void finishAnalysis(t_pbc *pbc);
 266         virtual void finishAnalysis(int nframes);
 267
 268         //! Routine to write output, that is additional over the built-in
 269         virtual void writeOutput();
 270
 271     private:
 272
 273         std::string                                                 fnNdx_;
 274
 275         std::vector< int >                                          index;
 276         std::vector< std::vector < RVec > >                         trajectory;
 277         Selection                                                   selec;
 278         int                                                         frames              = 0;
 279
 280         real                                                        **w_rls;
 281
 282         int                                                         Fitng_ePBC;
 283         // Copy and assign disallowed by base.
 284 };
 285
 286 Fitng::Fitng(): TrajectoryAnalysisModule()
 287 {
 288 }
 289
 290 Fitng::~Fitng()
 291 {
 292 }
 293
 294 void
 295 Fitng::initOptions(IOptionsContainer          *options,
 296                      TrajectoryAnalysisSettings *settings)
 297 {
 298     static const char *const desc[] = {
 299         "[THISMODULE] to be done"
 300     };
 301     // Add the descriptive text (program help text) to the options
 302     settings->setHelpText(desc);
 303     // Add option for selecting a subset of atoms
 304     options->addOption(SelectionOption("select")
 305                            .store(&selec).required()
 306                            .description("Atoms that are considered as part of the excluded volume"));
 307     // Add option for output file name
 308     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 309                             .store(&fnNdx_).defaultBasename("Fitng")
 310                             .description("Index file from the Fitng"));
 311     // Add option for etalon_frame constant
 312     //options->addOption(gmx::IntegerOption("dms")
 313     //                        .store(&domain_min_size)
 314     //                        .description("minimum domain size"));
 315     // Add option for epsi constant
 316     //options->addOption(DoubleOption("epsilon")
 317     //                        .store(&epsi)
 318     //                        .description("thermal vibrations' constant"));
 319     // Add option for delta constant
 320     //options->addOption(DoubleOption("delta")
 321     //                        .store(&delta)
 322     //                        .description("domain membership probability"));
 323     // Control input settings
 324     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 325     settings->setPBC(true);
 326 }
 327
 328 void
 329 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
 330                       const TopologyInformation        &top)
 331 {
 332     Fitng_ePBC = top.ePBC();
 333 }
 334
 335 void
 336 Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 337                              const t_trxframe                       &fr)
 338 {
 339     t_pbc  pbc;
 340     t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
 341     matrix boxx;
 342     copy_mat(fr.box, boxx);
 343     if (ppbc != NULL) {
 344         set_pbc(ppbc, Fitng_ePBC, boxx);
 345     }
 346     ConstArrayRef< int > atomind  = selec.atomIndices();
 347     index.resize(0);
 348     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 349         index.push_back(*ai);
 350     }
 351     trajectory.resize(1);
 352     trajectory.back().resize(selec.atomCount());
 353
 354     for (int i = 0; i < selec.atomCount(); i++) {
 355         trajectory.back()[i] = fr.x[index[i]];
 356     }
 357 }
 358
 359 void
 360 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 361                       TrajectoryAnalysisModuleData *pdata)
 362 {
 363     trajectory.resize(trajectory.size() + 1);
 364     trajectory.back().resize(index.size());
 365     for (int i = 0; i < index.size(); i++) {
 366         trajectory.back()[i] = fr.x[index[i]];
 367     }
 368 }
 369
 370 //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'
 371
 372 void
 373 Fitng::finishAnalysis(int nframes)
 374 {
 375     std::vector< std::pair< int, int > > pairs;
 376     for (int i = 0; i < index.size(); i++) {
 377         pairs.push_back(std::make_pair(i, i));
 378     }
 379
 380     double dist1 = 0, dist2 = 0;
 381
 382     for (int i = 0; i < index.size(); i++) {
 383         dist1 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory.back()[i][0], trajectory.back()[i][1], trajectory.back()[i][2], 0, 0, 0, 0, 0, 0);
 384     }
 385
 386     std::cout << "\n\n\n";
 387     std::cout << "\n\n\n" << "old dist = " << dist1 << "\n";
 388
 389     rvec *traj1, *traj2;
 390     snew(traj1, index.size());
 391     snew(traj2, index.size());
 392     for (int i = 0; i < index.size(); i++) {
 393         copy_rvec(trajectory[0][i].as_vec(), traj1[i]);
 394         copy_rvec(trajectory.back()[i].as_vec(), traj2[i]);
 395     }
 396
 397     //
 398     //  My Fit
 399     //
 400
 401     MyFitNew(trajectory[0], trajectory.back(), pairs);
 402     for (int i = 0; i < index.size(); i++) {
 403         dist2 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory.back()[i][0], trajectory.back()[i][1], trajectory.back()[i][2], 0, 0, 0, 0, 0, 0);
 404     }
 405     std::cout << "\n" << "my fit dist = " << dist2 << "\n";
 406
 407
 408     dist2 = 0;
 409     std::vector< real > w_rls;
 410     w_rls.resize(index.size(), 1);
 411     reset_x(index.size(), NULL, index.size(), NULL, traj1, &w_rls[0]);
 412     reset_x(index.size(), NULL, index.size(), NULL, traj2, &w_rls[0]);
 413     //
 414     //  Old Fit
 415     //
 416     do_fit(index.size(), &w_rls[0], traj1, traj2);
 417
 418     for (int i = 0; i < index.size(); i++) {
 419         dist2 += F(traj1[i][0], traj1[i][1], traj1[i][2], traj2[i][0], traj2[i][1], traj2[i][2], 0, 0, 0, 0, 0, 0);
 420     }
 421     std::cout << "\n" << "old fit dist = " << dist2 << "\n";
 422
 423     sfree(traj1);
 424     sfree(traj2);
 425 }
 426
 427 void
 428 Fitng::writeOutput()
 429 {
 430
 431 }
 432
 433
 434 /*! \brief
 435  * The main function for the analysis template.
 436  */
 437 int
 438 main(int argc, char *argv[])
 439 {
 440     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
 441 }