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 long double F (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double Fx (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double Fy (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double Fz (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double FA (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double FB (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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 long double FC (long double aix, long double aiy, long double aiz, long double bix, long double biy, long double biz, long double x, long double y, long double z, long double A, long double B, long 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, long double x, long double y, long double z, long double A, long double B, long 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     long double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0;
 167     long double f21 = 0, f22 = 0, f23 = 0, f24 = 0, f25 = 0, f26 = 0;
 168     long 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     long 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         if (f1 != 0) {
 186             while (true) {
 187                 f2 = 0;
 188                 for (int i = 0; i < pairs.size(); i++) {
 189                     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);
 190                     /*
 191                     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);
 192                     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);
 193                     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);
 194                     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);
 195                     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);
 196                     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);
 197                     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);
 198                     */
 199                 }
 200                 /*count++;
 201                 if (count % 100 == 0) {
 202                     std::cout << "   " << f1 << " " << f2 << "\n";
 203                 }*/
 204                 if (f2 < f1) {
 205                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
 206                     //x -= l1 * fx; y -= l2 * fy; z -= l3 * fz; A -= l4 * fa; B -= l5 * fb; C -= l6 * fc;
 207                     fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
 208                     break;
 209                 } else {
 210                     l *= 0.85;
 211                     /*
 212                     if (f21 > f1) {
 213                         l1 /= 2;
 214                     }
 215                     if (f22 > f1) {
 216                         l2 /= 2;
 217                     }
 218                     if (f23 > f1) {
 219                         l3 /= 2;
 220                     }
 221                     if (f24 > f1) {
 222                         l4 /= 2;
 223                     }
 224                     if (f25 > f1) {
 225                         l5 /= 2;
 226                     }
 227                     if (f26 > f1) {
 228                         l6 /= 2;
 229                     }
 230                     */
 231                 }
 232             }
 233         }
 234         if (f1 - f2 >= 0 && f1 - f2 < 0.00001) {
 235             break;
 236         }
 237         f1 = 0; f2 = 0;
 238     }
 239     ApplyFit(b, x, y, z, A, B, C);
 240 }
 241
 242 class Fitng : public TrajectoryAnalysisModule
 243 {
 244     public:
 245
 246         Fitng();
 247         virtual ~Fitng();
 248
 249         //! Set the options and setting
 250         virtual void initOptions(IOptionsContainer          *options,
 251                                  TrajectoryAnalysisSettings *settings);
 252
 253         //! First routine called by the analysis framework
 254         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 255         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 256                                   const TopologyInformation        &top);
 257
 258         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 259                                          const t_trxframe                   &fr);
 260
 261         //! Call for each frame of the trajectory
 262         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 263         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 264                                   TrajectoryAnalysisModuleData *pdata);
 265
 266         //! Last routine called by the analysis framework
 267         // virtual void finishAnalysis(t_pbc *pbc);
 268         virtual void finishAnalysis(int nframes);
 269
 270         //! Routine to write output, that is additional over the built-in
 271         virtual void writeOutput();
 272
 273     private:
 274
 275         std::string                                                 fnNdx_;
 276
 277         std::vector< int >                                          index;
 278         std::vector< std::vector < RVec > >                         trajectory;
 279         Selection                                                   selec;
 280         int                                                         frames              = 0;
 281
 282         real                                                        **w_rls;
 283
 284         int                                                         Fitng_ePBC;
 285         // Copy and assign disallowed by base.
 286 };
 287
 288 Fitng::Fitng(): TrajectoryAnalysisModule()
 289 {
 290 }
 291
 292 Fitng::~Fitng()
 293 {
 294 }
 295
 296 void
 297 Fitng::initOptions(IOptionsContainer          *options,
 298                      TrajectoryAnalysisSettings *settings)
 299 {
 300     static const char *const desc[] = {
 301         "[THISMODULE] to be done"
 302     };
 303     // Add the descriptive text (program help text) to the options
 304     settings->setHelpText(desc);
 305     // Add option for selecting a subset of atoms
 306     options->addOption(SelectionOption("select")
 307                            .store(&selec).required()
 308                            .description("Atoms that are considered as part of the excluded volume"));
 309     // Add option for output file name
 310     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 311                             .store(&fnNdx_).defaultBasename("Fitng")
 312                             .description("Index file from the Fitng"));
 313     // Add option for etalon_frame constant
 314     //options->addOption(gmx::IntegerOption("dms")
 315     //                        .store(&domain_min_size)
 316     //                        .description("minimum domain size"));
 317     // Add option for epsi constant
 318     //options->addOption(DoubleOption("epsilon")
 319     //                        .store(&epsi)
 320     //                        .description("thermal vibrations' constant"));
 321     // Add option for delta constant
 322     //options->addOption(DoubleOption("delta")
 323     //                        .store(&delta)
 324     //                        .description("domain membership probability"));
 325     // Control input settings
 326     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 327     settings->setPBC(true);
 328 }
 329
 330 void
 331 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
 332                       const TopologyInformation        &top)
 333 {
 334     Fitng_ePBC = top.ePBC();
 335 }
 336
 337 void
 338 Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 339                              const t_trxframe                       &fr)
 340 {
 341     t_pbc  pbc;
 342     t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
 343     matrix boxx;
 344     copy_mat(fr.box, boxx);
 345     if (ppbc != NULL) {
 346         set_pbc(ppbc, Fitng_ePBC, boxx);
 347     }
 348     ConstArrayRef< int > atomind  = selec.atomIndices();
 349     index.resize(0);
 350     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 351         index.push_back(*ai);
 352     }
 353     trajectory.resize(1);
 354     trajectory.back().resize(selec.atomCount());
 355
 356     for (int i = 0; i < selec.atomCount(); i++) {
 357         trajectory.back()[i] = fr.x[index[i]];
 358     }
 359 }
 360
 361 void
 362 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 363                       TrajectoryAnalysisModuleData *pdata)
 364 {
 365     trajectory.resize(trajectory.size() + 1);
 366     trajectory.back().resize(index.size());
 367     for (int i = 0; i < index.size(); i++) {
 368         trajectory.back()[i] = fr.x[index[i]];
 369     }
 370 }
 371
 372 // 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'
 373 // Fitng -s '/home/toluk/Develop/FitSamples/cube_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube.pdb' -select 'name NA'
 374 // Fitng -s '/home/toluk/Develop/FitSamples/cube_with_termal_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_termal.pdb' -select 'name NA'
 375 // 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'
 376 // 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'
 377 // Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
 378 // Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_2nd_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
 379 // Fitng -s '/home/toluk/Develop/FitSamples/rubber_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/rubber.pdb' -select 'name NA'
 380 // rubber_1st_fr
 381
 382 void
 383 Fitng::finishAnalysis(int nframes)
 384 {
 385     std::vector< std::pair< int, int > > pairs;
 386     for (int i = 0; i < index.size(); i++) {
 387         pairs.push_back(std::make_pair(i, i));
 388     }
 389     double dist1 = 0, dist2 = 0;
 390     std::vector< real > w_rls;
 391     w_rls.resize(index.size(), 1);
 392     rvec *traj1, **traj2;
 393     snew(traj1, index.size());
 394     snew(traj2, trajectory.size());
 395     for (int i = 0; i < index.size(); i++) {
 396         copy_rvec(trajectory[0][i].as_vec(), traj1[i]);
 397     }
 398     reset_x(index.size(), NULL, index.size(), NULL, traj1, &w_rls[0]);
 399     for (int j = 0; j < trajectory.size(); j++) {
 400         snew(traj2[j], index.size());
 401         for (int i = 0; i < index.size(); i++) {
 402             copy_rvec(trajectory[j][i].as_vec(), traj2[j][i]);
 403         }
 404     }
 405
 406     for (int j = 1; j < trajectory.size(); j++) {
 407
 408         for (int i = 0; i < index.size(); i++) {
 409             dist1 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory[j][i][0], trajectory[j][i][1], trajectory[j][i][2], 0, 0, 0, 0, 0, 0);
 410             //std::cout << F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory[j][i][0], trajectory[j][i][1], trajectory[j][i][2], 0, 0, 0, 0, 0, 0);
 411         }
 412         std::cout << "\nbasic dist = " << dist1 / index.size() << "\n";
 413
 414         //
 415         //  My Fit
 416         //
 417
 418         MyFitNew(trajectory[0], trajectory[j], pairs);
 419         for (int i = 0; i < index.size(); i++) {
 420             dist2 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory[j][i][0], trajectory[j][i][1], trajectory[j][i][2], 0, 0, 0, 0, 0, 0);
 421         }
 422         std::cout << "my fit dist = " << dist2 / index.size()<< "\n";
 423
 424         dist2 = 0;
 425
 426         reset_x(index.size(), NULL, index.size(), NULL, traj2[j], &w_rls[0]);
 427         //
 428         //  Old Fit
 429         //
 430
 431         /*for (int i = 0; i < index.size(); i++) {
 432             std::cout << traj2[j][i][0] << " " << traj2[j][i][1] << " " << traj2[j][i][2] << "\n";
 433         }*/
 434         do_fit(index.size(), &w_rls[0], traj1, traj2[j]);
 435         /*for (int i = 0; i < index.size(); i++) {
 436             std::cout << traj2[j][i][0] << " " << traj2[j][i][1] << " " << traj2[j][i][2] << "\n";
 437         }*/
 438
 439         for (int i = 0; i < index.size(); i++) {
 440             dist2 += F(traj1[i][0], traj1[i][1], traj1[i][2], traj2[j][i][0], traj2[j][i][1], traj2[j][i][2], 0, 0, 0, 0, 0, 0);
 441             //std::cout << F(traj1[i][0], traj1[i][1], traj1[i][2], traj2[j][i][0], traj2[j][i][1], traj2[j][i][2], 0, 0, 0, 0, 0, 0);
 442         }
 443         std::cout << "old fit dist = " << dist2 / index.size()<< "\n";
 444
 445         dist1 = 0;
 446         dist2 = 0;
 447
 448     }
 449
 450     int co = 0;
 451     freopen ("/home/toluk/Develop/FitSamples/old_fit_result.pdb","w",stdout);
 452     for (int l = 1; l < trajectory.size(); l++) {
 453         for (int i = 0; i < index.size(); i++) {
 454             //cout << D*i + (float)l * 6;
 455             printf("ATOM  %5d  NA   NA %5d    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s  \n", co, co, (traj2[l][i][0] * 10), (traj2[l][i][1] * 10), (traj2[l][i][2] * 10), 1.0, 20.0, "    ", " ");
 456             co++;
 457         }
 458         co = 0;
 459         printf("ENDMDL\n");
 460     }
 461     printf("END\n");
 462     freopen ("/home/toluk/Develop/FitSamples/new_fit_result.pdb","w",stdout);
 463     for (int l = 1; l < trajectory.size(); l++) {
 464         for (int i = 0; i < index.size(); i++) {
 465             //cout << D*i + (float)l * 6;
 466             printf("ATOM  %5d  NA   NA %5d    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s  \n", co, co, (trajectory[l][i][0] * 10), (trajectory[l][i][1] * 10), (trajectory[l][i][2] * 10), 1.0, 20.0, "    ", " ");
 467             co++;
 468         }
 469         co = 0;
 470         printf("ENDMDL\n");
 471     }
 472     printf("END\n");
 473     fclose(stdout);
 474
 475     sfree(traj1);
 476     for (int i = 0; i < trajectory.size(); i++) {
 477         sfree(traj2[i]);
 478     }
 479     sfree(traj2);
 480 }
 481
 482 void
 483 Fitng::writeOutput()
 484 {
 485
 486 }
 487
 488
 489 /*! \brief
 490  * The main function for the analysis template.
 491  */
 492 int
 493 main(int argc, char *argv[])
 494 {
 495     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
 496 }