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             std::cout << f1 << " " << f2 << "\n";
 200             if (f2 < f1) {
 201                 x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 202                 //x -= l1 * fx; y -= l2 * fy; z -= l3 * fz; A -= l4 * fa; B -= l5 * fb; C -= l6 * fc;
 203                 fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 204                 break;
 205             } else {
 206                 l *= 0.85;
 207                 /*
 208                 if (f21 > f1) {
 209                     l1 /= 2;
 210                 }
 211                 if (f22 > f1) {
 212                     l2 /= 2;
 213                 }
 214                 if (f23 > f1) {
 215                     l3 /= 2;
 216                 }
 217                 if (f24 > f1) {
 218                     l4 /= 2;
 219                 }
 220                 if (f25 > f1) {
 221                     l5 /= 2;
 222                 }
 223                 if (f26 > f1) {
 224                     l6 /= 2;
 225                 }
 226                 */
 227             }
 228         }
 229         if (f1 - f2 > 0 && f1 - f2 < 0.00001) {
 230             break;
 231         }
 232         f1 = 0; f2 = 0;
 233     }
 234     ApplyFit(b, x, y, z, A, B, C);
 235 }
 236
 237 class Fitng : public TrajectoryAnalysisModule
 238 {
 239     public:
 240
 241         Fitng();
 242         virtual ~Fitng();
 243
 244         //! Set the options and setting
 245         virtual void initOptions(IOptionsContainer          *options,
 246                                  TrajectoryAnalysisSettings *settings);
 247
 248         //! First routine called by the analysis framework
 249         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 250         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 251                                   const TopologyInformation        &top);
 252
 253         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 254                                          const t_trxframe                   &fr);
 255
 256         //! Call for each frame of the trajectory
 257         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 258         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 259                                   TrajectoryAnalysisModuleData *pdata);
 260
 261         //! Last routine called by the analysis framework
 262         // virtual void finishAnalysis(t_pbc *pbc);
 263         virtual void finishAnalysis(int nframes);
 264
 265         //! Routine to write output, that is additional over the built-in
 266         virtual void writeOutput();
 267
 268     private:
 269
 270         std::string                                                 fnNdx_;
 271
 272         std::vector< int >                                          index;
 273         std::vector< std::vector < RVec > >                         trajectory;
 274         Selection                                                   selec;
 275         int                                                         frames              = 0;
 276
 277         real                                                        **w_rls;
 278
 279         int                                                         Fitng_ePBC;
 280         // Copy and assign disallowed by base.
 281 };
 282
 283 Fitng::Fitng(): TrajectoryAnalysisModule()
 284 {
 285 }
 286
 287 Fitng::~Fitng()
 288 {
 289 }
 290
 291 void
 292 Fitng::initOptions(IOptionsContainer          *options,
 293                      TrajectoryAnalysisSettings *settings)
 294 {
 295     static const char *const desc[] = {
 296         "[THISMODULE] to be done"
 297     };
 298     // Add the descriptive text (program help text) to the options
 299     settings->setHelpText(desc);
 300     // Add option for selecting a subset of atoms
 301     options->addOption(SelectionOption("select")
 302                            .store(&selec).required()
 303                            .description("Atoms that are considered as part of the excluded volume"));
 304     // Add option for output file name
 305     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 306                             .store(&fnNdx_).defaultBasename("Fitng")
 307                             .description("Index file from the Fitng"));
 308     // Add option for etalon_frame constant
 309     //options->addOption(gmx::IntegerOption("dms")
 310     //                        .store(&domain_min_size)
 311     //                        .description("minimum domain size"));
 312     // Add option for epsi constant
 313     //options->addOption(DoubleOption("epsilon")
 314     //                        .store(&epsi)
 315     //                        .description("thermal vibrations' constant"));
 316     // Add option for delta constant
 317     //options->addOption(DoubleOption("delta")
 318     //                        .store(&delta)
 319     //                        .description("domain membership probability"));
 320     // Control input settings
 321     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 322     settings->setPBC(true);
 323 }
 324
 325 void
 326 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
 327                       const TopologyInformation        &top)
 328 {
 329     Fitng_ePBC = top.ePBC();
 330 }
 331
 332 void
 333 Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 334                              const t_trxframe                       &fr)
 335 {
 336     t_pbc  pbc;
 337     t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
 338     matrix boxx;
 339     copy_mat(fr.box, boxx);
 340     if (ppbc != NULL) {
 341         set_pbc(ppbc, Fitng_ePBC, boxx);
 342     }
 343     ConstArrayRef< int > atomind  = selec.atomIndices();
 344     index.resize(0);
 345     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 346         index.push_back(*ai);
 347     }
 348     trajectory.resize(1);
 349     trajectory.back().resize(selec.atomCount());
 350
 351     for (int i = 0; i < selec.atomCount(); i++) {
 352         trajectory.back()[i] = fr.x[index[i]];
 353     }
 354 }
 355
 356 void
 357 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 358                       TrajectoryAnalysisModuleData *pdata)
 359 {
 360     trajectory.resize(trajectory.size() + 1);
 361     trajectory.back().resize(index.size());
 362     for (int i = 0; i < index.size(); i++) {
 363         trajectory.back()[i] = fr.x[index[i]];
 364     }
 365 }
 366
 367 //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'
 368
 369 void
 370 Fitng::finishAnalysis(int nframes)
 371 {
 372     std::vector< std::pair< int, int > > pairs;
 373     for (int i = 0; i < index.size(); i++) {
 374         pairs.push_back(std::make_pair(i, i));
 375     }
 376
 377     double dist1 = 0, dist2 = 0;
 378
 379     for (int i = 0; i < index.size(); i++) {
 380         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);
 381     }
 382
 383     std::cout << "\n\n\n";
 384
 385     MyFitNew(trajectory[0], trajectory.back(), pairs);
 386
 387     for (int i = 0; i < index.size(); i++) {
 388         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);
 389     }
 390
 391     std::cout << "\n\n\n" << "old dist = " << dist1 << "\nnew dist = " << dist2 << "\n\n";
 392
 393 }
 394
 395 void
 396 Fitng::writeOutput()
 397 {
 398
 399 }
 400
 401
 402 /*! \brief
 403  * The main function for the analysis template.
 404  */
 405 int
 406 main(int argc, char *argv[])
 407 {
 408     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
 409 }