src/spirals.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2011,2012,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 #include <omp.h>
  36 #include <iostream>
  37 #include <math.h>
  38
  39 #include <gromacs/trajectoryanalysis.h>
  40 #include <gromacs/math/do_fit.h>
  41 #include <gromacs/utility/smalloc.h>
  42 #include "gromacs/selection/selection.h"
  43 #include "gromacs/selection/selectionoption.h"
  44
  45 using namespace gmx;
  46
  47 struct kernel_maxima {
  48     RVec x;
  49     RVec p;
  50     std::vector< RVec > krnl;
  51 };
  52
  53 long double Fx (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
  54     long double ret = 0;
  55     for (int i = 0; i < x.size(); i++) {
  56         ret +=
  57            sqrt (   pow (p2 * (x[i][2] - z0) - p3 * (x[i][1] - y0), 2) +
  58                         pow (p3 * (x[i][0] - x0) - p1 * (x[i][2] - z0), 2) +
  59                         pow (p1 * (x[i][1] - y0) - p2 * (x[i][0] - x0), 2)) /
  60            sqrt (p1 * p1 + p2 * p2 + p3 * p3);
  61     }
  62     return ret;
  63 }
  64
  65 long double fx0 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
  66     long double ret = 0;
  67     for (int i = 0; i < x.size(); i++) {
  68         ret +=
  69            (2 * p2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) + 2 * p3 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]))) /
  70            (2 *     sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  71                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  72                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
  73                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
  74     }
  75     return ret;
  76 }
  77
  78 long double fy0 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
  79     long double ret = 0;
  80     for (int i = 0; i < x.size(); i++) {
  81         ret +=
  82            -(2 * p1 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) - 2 * p3 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
  83             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  84                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  85                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
  86                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
  87     }
  88     return ret;
  89 }
  90
  91 long double fz0 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
  92     long double ret = 0;
  93     for (int i = 0; i < x.size(); i++) {
  94         ret +=
  95            -(2 * p1 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) + 2 * p2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
  96             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  97                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  98                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
  99                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
 100     }
 101     return ret;
 102 }
 103
 104 long double fp1 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
 105     long double ret = 0;
 106     for (int i = 0; i < x.size(); i++) {
 107         ret +=
 108            -(2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) * (y0 - x[i][1]) + 2 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) * (z0 - x[i][2])) /
 109             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 110                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 111                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 112                     sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 113             (p1 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 114                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 115                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 116                     pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 117     }
 118     return ret;
 119 }
 120
 121 long double fp2 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
 122     long double ret = 0;
 123     for (int i = 0; i < x.size(); i++) {
 124         ret +=
 125            (2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) * (x0 - x[i][0]) - 2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2])) * (z0 - x[i][2])) /
 126            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 127                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 128                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 129                    sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 130            (p2 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 131                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 132                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 133                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 134     }
 135     return ret;
 136 }
 137
 138 long double fp3 (long double x0, long double y0, long double z0, long double p1, long double p2, long double p3, std::vector< RVec > x) {
 139     long double ret = 0;
 140     for (int i = 0; i < x.size(); i++) {
 141         ret +=
 142            (2 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) * (x0 - x[i][0]) + 2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2])) * (y0 - x[i][1])) /
 143            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 144                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 145                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 146                    sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 147            (p3 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 148                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 149                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 150                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 151     }
 152     return ret;
 153 }
 154
 155 void linear_kernel_search (long double &x0, long double &y0, long double &z0, long double &p1, long double &p2, long double &p3, std::vector< RVec > x, long double epsi) {
 156     long double FX = 0, FX0 = 0, FY0 = 0, FZ0 = 0, FP1 = 0, FP2 = 0, FP3 = 0;
 157     long double L1, L2, L3, L4, L5, L6;
 158     //int count = 0;
 159     while (true) {
 160         FX = Fx(x0, y0, z0, p1, p2, p3, x);
 161         FX0 = fx0(x0, y0, z0, p1, p2, p3, x);
 162         FY0 = fy0(x0, y0, z0, p1, p2, p3, x);
 163         FZ0 = fz0(x0, y0, z0, p1, p2, p3, x);
 164         FP1 = fp1(x0, y0, z0, p1, p2, p3, x);
 165         FP2 = fp2(x0, y0, z0, p1, p2, p3, x);
 166         FP3 = fp3(x0, y0, z0, p1, p2, p3, x);
 167
 168         L1 = 1;
 169         while (Fx(x0 - L1 * FX0, y0, z0, p1, p2, p3, x) - FX > 0) {
 170             L1 /= 2;
 171             if (x0 - L1 * FX0 < epsi) {
 172                 L1 = 0;
 173             }
 174         }
 175         std::cout << FX - Fx(x0 - L1 * FX0, y0, z0, p1, p2, p3, x) << " ";
 176         L2 = 1;
 177         while (Fx(x0, y0 - L2 * FY0, z0, p1, p2, p3, x) - FX > 0) {
 178             L2 /= 2;
 179             if (y0 - L2 * FY0 < epsi) {
 180                 L2 = 0;
 181             }
 182         }
 183         std::cout << FX - Fx(x0, y0 - L2 * FY0, z0, p1, p2, p3, x) << " ";
 184         L3 = 1;
 185         while (Fx(x0, y0, z0 - L3 * FZ0, p1, p2, p3, x) - FX > 0) {
 186             L3 /= 2;
 187             if (z0 - L3 * FZ0 < epsi) {
 188                 L3 = 0;
 189             }
 190         }
 191         std::cout << FX - Fx(x0, y0, z0 - L3 * FZ0, p1, p2, p3, x) << " ";
 192         L4 = 1;
 193         while (Fx(x0, y0, z0, p1 - L4 * FP1, p2, p3, x) - FX > 0) {
 194             L4 /= 2;
 195             if (p1 - L4 * FP1 < epsi) {
 196                 L4 = 0;
 197             }
 198         }
 199         std::cout << FX - Fx(x0, y0, z0, p1 - L4 * FP1, p2, p3, x) << " ";
 200         L5 = 1;
 201         while (Fx(x0, y0, z0, p1, p2 - L5 * FP2, p3, x) - FX > 0) {
 202             L5 /= 2;
 203             if (p2 - L5 * FP2 < epsi) {
 204                 L5 = 0;
 205             }
 206         }
 207         std::cout << FX - Fx(x0, y0, z0, p1, p2 - L5 * FP2, p3, x) << " ";
 208         L6 = 1;
 209         while (Fx(x0, y0, z0, p1, p2, p3 - L6 * FP3, x) - FX > 0) {
 210             L6 /= 2;
 211             if (p3 - L6 * FP3 < epsi) {
 212                 L6 = 0;
 213             }
 214         }
 215         std::cout << FX - Fx(x0, y0, z0, p1, p2, p3 - L6 * FP3, x) << " ";
 216         std::cout << FX - Fx(x0 - L1 * FX0, y0 - L2 * FY0, z0 - L3 * FZ0, p1 - L4 * FP1, p2 - L5 * FP2, p3 - L6 * FP3, x) << "\n";
 217         if (FX - Fx(x0 - L1 * FX0, y0 - L2 * FY0, z0 - L3 * FZ0, p1 - L4 * FP1, p2 - L5 * FP2, p3 - L6 * FP3, x) > epsi) {
 218             x0 -= L1 * FX0;
 219             y0 -= L2 * FY0;
 220             z0 -= L3 * FZ0;
 221             p1 -= L4 * FP1;
 222             p2 -= L5 * FP2;
 223             p3 -= L6 * FP3;
 224         } else {
 225             break;
 226         }
 227         //count++;
 228     }
 229     std::cout << "\n";
 230     //std::cout << count << "\n";
 231 }
 232
 233 RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p3, RVec x) {
 234     double lambda = (p1 * (x[0] - x0) + p2 * (x[1] - y0) + p3 * (x[2] - z0)) / (p1 * p1 + p2 * p2 + p3 * p3);
 235     RVec pro;
 236     pro[0] = x0 + p1 * lambda;
 237     pro[1] = y0 + p2 * lambda;
 238     pro[2] = z0 + p3 * lambda;
 239     return pro;
 240 }
 241
 242 double left_right_turn (RVec a, RVec b, RVec c) {
 243     return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + b[0] * c[1] * a[2] -
 244                 (a[2] * b[1] * c[0] + a[0] * b[2] * c[1] + a[1] * b[0] * c[2]);
 245 }
 246
 247 /*! \brief
 248  * Template class to serve as a basis for user analysis tools.
 249  */
 250 class Spirals : public TrajectoryAnalysisModule
 251 {
 252     public:
 253         Spirals();
 254
 255         virtual void initOptions(IOptionsContainer          *options,
 256                                  TrajectoryAnalysisSettings *settings);
 257         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 258                                   const TopologyInformation        &top);
 259
 260         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 261                                   TrajectoryAnalysisModuleData *pdata);
 262
 263         virtual void finishAnalysis(int nframes);
 264         virtual void writeOutput();
 265
 266     private:
 267         class ModuleData;
 268
 269         std::string                      fnNdx_;
 270         double                           cutoff_;
 271         Selection                        refsel_;
 272
 273         AnalysisNeighborhood             nb_;
 274
 275         AnalysisData                     data_;
 276         AnalysisDataAverageModulePointer avem_;
 277
 278         SelectionList                           sel_;
 279         int                                     frames      = 0;
 280         double                                  epsi        = 0.00001;
 281
 282         std::vector< std::vector< RVec > >                  monomers;
 283         std::vector< kernel_maxima >                        kernel;
 284         std::vector< std::vector< std::vector< int > > >    circles;
 285 };
 286
 287 Spirals::Spirals()
 288     : cutoff_(0.0)
 289 {
 290     registerAnalysisDataset(&data_, "avedist");
 291 }
 292
 293 void
 294 Spirals::initOptions(IOptionsContainer          *options,
 295                               TrajectoryAnalysisSettings *settings)
 296 {
 297     static const char *const desc[] = {
 298         "Analysis tool for finding molecular core."
 299     };
 300
 301     // Add the descriptive text (program help text) to the options
 302     settings->setHelpText(desc);
 303     // Add option for output file name
 304     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 305                             .store(&fnNdx_).defaultBasename("rcore")
 306                             .description("Index file from the rcore"));
 307     // Add option for selection list
 308     options->addOption(SelectionOption("select").storeVector(&sel_)
 309                            .required().dynamicMask().multiValue()
 310                            .description("Position to calculate distances for"));
 311 }
 312
 313 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test.ndx' -on '/home/toluk/Develop/samples/reca_rd/core.ndx'
 314
 315 void
 316 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 317                                const TopologyInformation         & /*top*/)
 318 {
 319     kernel.resize(0);
 320     circles.resize(0);
 321     monomers.resize(0);
 322 }
 323
 324 void
 325 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 326                                TrajectoryAnalysisModuleData *pdata)
 327 {
 328     const SelectionList &sel = pdata->parallelSelections(sel_);
 329     std::vector< RVec > temp;
 330     temp.resize(sel_.size());
 331     for (int i = 0; i < sel.size(); i++) {
 332         copy_rvec(sel[i].position(0).x(), temp[i]);
 333     }
 334
 335     monomers.resize(monomers.size() + 1);
 336     for (int i = 0; i < sel.size(); i++) {
 337         monomers.back().push_back(temp[i]);
 338     }
 339
 340     RVec mid, arrow;
 341
 342     mid[0] = 0;
 343     mid[1] = 0;
 344     mid[2] = 0;
 345
 346     arrow[0] = 0;
 347     arrow[1] = 0;
 348     arrow[2] = 0;
 349
 350     for (int i = 0; i < sel.size(); i++) {
 351         rvec_inc(temp[i], mid);
 352     }
 353     mid[0] /= sel.size();
 354     mid[1] /= sel.size();
 355     mid[2] /= sel.size();
 356     rvec_sub(temp.back(), temp.front(), arrow);
 357
 358     long double t1, t2, t3, t4, t5, t6;
 359     t1 = mid[0];
 360     t2 = mid[1];
 361     t3 = mid[2];
 362     t4 = arrow[0];
 363     t5 = arrow[1];
 364     t6 = arrow[2];
 365
 366     //linear_kernel_search(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp, epsi); //изменить формат функции для изменения значений переменных
 367     linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
 368
 369     //std::cout << t1 << " " << t2 << " " << t3 << " " << t4 << " " << t5 << " " << t6 << "\n";
 370
 371     mid[0] = t1;
 372     mid[1] = t2;
 373     mid[2] = t3;
 374     arrow[0] = t4;
 375     arrow[1] = t5;
 376     arrow[2] = t6;
 377
 378     kernel.resize(kernel.size() + 1);
 379     kernel.back().x = mid;
 380     kernel.back().p = arrow;
 381     for (int i = 0; i < sel.size(); i++) {
 382         kernel.back().krnl.push_back(kernel_pro(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp[i]));
 383     }
 384
 385     /*circles.resize(circles.size() + 1);
 386
 387     int prev_sign = -9999, change = 0;
 388     RVec a, b, c;
 389     rvec_sub(temp[0], kernel.back().krnl.front(), a);
 390     rvec_sub(kernel.back().krnl.front(), kernel.back().krnl.back(), b);
 391     for (int i = 1; i < sel.size(); i++) {
 392         rvec_sub(temp[i], kernel.back().krnl[i], c);
 393         if (std::signbit(left_right_turn(a, b, c)) != prev_sign) {
 394             change++;
 395             prev_sign == std::signbit(left_right_turn(a, b, c));
 396         }
 397         if (change != 3) {
 398             if (circles.back().size() == 0) {
 399                 circles.back().resize(1);
 400             }
 401             circles.back().back().push_back(i);
 402         } else {
 403             circles.back().resize(circles.back().size() + 1);
 404             circles.back().back().push_back(i);
 405             change = 1;
 406         }
 407     }*/
 408
 409     frames++;
 410 }
 411
 412 void
 413 Spirals::finishAnalysis(int /*nframes*/)
 414 {
 415
 416
 417
 418     FILE *file;
 419     file = std::fopen("linear_kernel.txt", "w+");
 420     for (int i = 0; i < kernel.size(); i++) {
 421         for (int j = 0; j < monomers[i].size(); j++) {
 422             std::fprintf(file, "%3.2f %3.2f %3.2f\n", monomers[i][j][0], monomers[i][j][1], monomers[i][j][2]);
 423         }
 424         for (int j = 0; j < kernel[i].krnl.size(); j++) {
 425             std::fprintf(file, "%3.2f %3.2f %3.2f\n", kernel[i].krnl[j][0], kernel[i].krnl[j][1], kernel[i].krnl[j][2]);
 426         }
 427         std::fprintf(file, "\n");
 428     }
 429     std::fclose(file);
 430
 431
 432
 433     /*
 434     FILE *file;
 435     file = std::fopen("spiral_dist.txt", "w+");
 436     for (int i = 0; i < spiral_dist.size(); i++) {
 437         for (int j = 0; j < spiral_dist[i].size(); j++) {
 438             std::fprintf(file, "%3.2Lf\n", spiral_dist[i][j]);
 439         }
 440         std::fprintf(file, "\n");
 441     }
 442     std::fclose(file);*/
 443 }
 444
 445 void
 446 Spirals::writeOutput()
 447 {
 448
 449 }
 450
 451 /*! \brief
 452  * The main function for the analysis template.
 453  */
 454 int
 455 main(int argc, char *argv[])
 456 {
 457     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 458 }