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         L2 = 1;
 176         while (Fx(x0, y0 - L2 * FY0, z0, p1, p2, p3, x) - FX > 0) {
 177             L2 /= 2;
 178             if (y0 - L2 * FY0 < epsi) {
 179                 L2 = 0;
 180             }
 181         }
 182         L3 = 1;
 183         while (Fx(x0, y0, z0 - L3 * FZ0, p1, p2, p3, x) - FX > 0) {
 184             L3 /= 2;
 185             if (z0 - L3 * FZ0 < epsi) {
 186                 L3 = 0;
 187             }
 188         }
 189         L4 = 1;
 190         while (Fx(x0, y0, z0, p1 - L4 * FP1, p2, p3, x) - FX > 0) {
 191             L4 /= 2;
 192             if (p1 - L4 * FP1 < epsi) {
 193                 L4 = 0;
 194             }
 195         }
 196         L5 = 1;
 197         while (Fx(x0, y0, z0, p1, p2 - L5 * FP2, p3, x) - FX > 0) {
 198             L5 /= 2;
 199             if (p2 - L5 * FP2 < epsi) {
 200                 L5 = 0;
 201             }
 202         }
 203         L6 = 1;
 204         while (Fx(x0, y0, z0, p1, p2, p3 - L6 * FP3, x) - FX > 0) {
 205             L6 /= 2;
 206             if (p3 - L6 * FP3 < epsi) {
 207                 L6 = 0;
 208             }
 209         }
 210         std::cout << FX - Fx(x0 - L1 * FX0, y0 - L2 * FY0, z0 - L3 * FZ0, p1 - L4 * FP1, p2 - L5 * FP2, p3 - L6 * FP3, x) << " ";
 211         if (FX - Fx(x0 - L1 * FX0, y0 - L2 * FY0, z0 - L3 * FZ0, p1 - L4 * FP1, p2 - L5 * FP2, p3 - L6 * FP3, x) > epsi) {
 212             x0 -= L1 * FX0;
 213             y0 -= L2 * FY0;
 214             z0 -= L3 * FZ0;
 215             p1 -= L4 * FP1;
 216             p2 -= L5 * FP2;
 217             p3 -= L6 * FP3;
 218         } else {
 219             break;
 220         }
 221         //count++;
 222     }
 223     std::cout << "\n";
 224     //std::cout << count << "\n";
 225 }
 226
 227 RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p3, RVec x) {
 228     double lambda = (p1 * (x[0] - x0) + p2 * (x[1] - y0) + p3 * (x[2] - z0)) / (p1 * p1 + p2 * p2 + p3 * p3);
 229     RVec pro;
 230     pro[0] = x0 + p1 * lambda;
 231     pro[1] = y0 + p2 * lambda;
 232     pro[2] = z0 + p3 * lambda;
 233     return pro;
 234 }
 235
 236 double left_right_turn (RVec a, RVec b, RVec c) {
 237     return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + b[0] * c[1] * a[2] -
 238                 (a[2] * b[1] * c[0] + a[0] * b[2] * c[1] + a[1] * b[0] * c[2]);
 239 }
 240
 241 /*! \brief
 242  * Template class to serve as a basis for user analysis tools.
 243  */
 244 class Spirals : public TrajectoryAnalysisModule
 245 {
 246     public:
 247         Spirals();
 248
 249         virtual void initOptions(IOptionsContainer          *options,
 250                                  TrajectoryAnalysisSettings *settings);
 251         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 252                                   const TopologyInformation        &top);
 253
 254         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 255                                   TrajectoryAnalysisModuleData *pdata);
 256
 257         virtual void finishAnalysis(int nframes);
 258         virtual void writeOutput();
 259
 260     private:
 261         class ModuleData;
 262
 263         std::string                      fnNdx_;
 264         double                           cutoff_;
 265         Selection                        refsel_;
 266
 267         AnalysisNeighborhood             nb_;
 268
 269         AnalysisData                     data_;
 270         AnalysisDataAverageModulePointer avem_;
 271
 272         SelectionList                           sel_;
 273         int                                     frames      = 0;
 274         double                                  epsi        = 0.00001;
 275
 276         std::vector< std::vector< RVec > >                  monomers;
 277         std::vector< kernel_maxima >                        kernel;
 278         std::vector< std::vector< std::vector< int > > >    circles;
 279 };
 280
 281 Spirals::Spirals()
 282     : cutoff_(0.0)
 283 {
 284     registerAnalysisDataset(&data_, "avedist");
 285 }
 286
 287 void
 288 Spirals::initOptions(IOptionsContainer          *options,
 289                               TrajectoryAnalysisSettings *settings)
 290 {
 291     static const char *const desc[] = {
 292         "Analysis tool for finding molecular core."
 293     };
 294
 295     // Add the descriptive text (program help text) to the options
 296     settings->setHelpText(desc);
 297     // Add option for output file name
 298     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 299                             .store(&fnNdx_).defaultBasename("rcore")
 300                             .description("Index file from the rcore"));
 301     // Add option for selection list
 302     options->addOption(SelectionOption("select").storeVector(&sel_)
 303                            .required().dynamicMask().multiValue()
 304                            .description("Position to calculate distances for"));
 305 }
 306
 307 // -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'
 308
 309 void
 310 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 311                                const TopologyInformation         & /*top*/)
 312 {
 313     kernel.resize(0);
 314     circles.resize(0);
 315     monomers.resize(0);
 316 }
 317
 318 void
 319 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 320                                TrajectoryAnalysisModuleData *pdata)
 321 {
 322     const SelectionList &sel = pdata->parallelSelections(sel_);
 323     std::vector< RVec > temp;
 324     temp.resize(sel_.size());
 325     for (int i = 0; i < sel.size(); i++) {
 326         copy_rvec(sel[i].position(0).x(), temp[i]);
 327     }
 328
 329     monomers.resize(monomers.size() + 1);
 330     for (int i = 0; i < sel.size(); i++) {
 331         monomers.back().push_back(temp[i]);
 332     }
 333
 334     RVec mid, arrow;
 335
 336     mid[0] = 0;
 337     mid[1] = 0;
 338     mid[2] = 0;
 339
 340     arrow[0] = 0;
 341     arrow[1] = 0;
 342     arrow[2] = 0;
 343
 344     for (int i = 0; i < sel.size(); i++) {
 345         rvec_inc(temp[i], mid);
 346     }
 347     mid[0] /= sel.size();
 348     mid[1] /= sel.size();
 349     mid[2] /= sel.size();
 350     rvec_sub(temp.back(), temp.front(), arrow);
 351
 352     long double t1, t2, t3, t4, t5, t6;
 353     t1 = mid[0];
 354     t2 = mid[1];
 355     t3 = mid[2];
 356     t4 = arrow[0];
 357     t5 = arrow[1];
 358     t6 = arrow[2];
 359
 360     //linear_kernel_search(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp, epsi); //изменить формат функции для изменения значений переменных
 361     linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
 362
 363     //std::cout << t1 << " " << t2 << " " << t3 << " " << t4 << " " << t5 << " " << t6 << "\n";
 364
 365     mid[0] = t1;
 366     mid[1] = t2;
 367     mid[2] = t3;
 368     arrow[0] = t4;
 369     arrow[1] = t5;
 370     arrow[2] = t6;
 371
 372     kernel.resize(kernel.size() + 1);
 373     kernel.back().x = mid;
 374     kernel.back().p = arrow;
 375     for (int i = 0; i < sel.size(); i++) {
 376         kernel.back().krnl.push_back(kernel_pro(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp[i]));
 377     }
 378
 379     /*circles.resize(circles.size() + 1);
 380
 381     int prev_sign = -9999, change = 0;
 382     RVec a, b, c;
 383     rvec_sub(temp[0], kernel.back().krnl.front(), a);
 384     rvec_sub(kernel.back().krnl.front(), kernel.back().krnl.back(), b);
 385     for (int i = 1; i < sel.size(); i++) {
 386         rvec_sub(temp[i], kernel.back().krnl[i], c);
 387         if (std::signbit(left_right_turn(a, b, c)) != prev_sign) {
 388             change++;
 389             prev_sign == std::signbit(left_right_turn(a, b, c));
 390         }
 391         if (change != 3) {
 392             if (circles.back().size() == 0) {
 393                 circles.back().resize(1);
 394             }
 395             circles.back().back().push_back(i);
 396         } else {
 397             circles.back().resize(circles.back().size() + 1);
 398             circles.back().back().push_back(i);
 399             change = 1;
 400         }
 401     }*/
 402
 403     frames++;
 404 }
 405
 406 void
 407 Spirals::finishAnalysis(int /*nframes*/)
 408 {
 409
 410
 411
 412     FILE *file;
 413     file = std::fopen("linear_kernel.txt", "w+");
 414     for (int i = 0; i < kernel.size(); i++) {
 415         for (int j = 0; j < monomers[i].size(); j++) {
 416             std::fprintf(file, "%3.2f %3.2f %3.2f\n", monomers[i][j][0], monomers[i][j][1], monomers[i][j][2]);
 417         }
 418         for (int j = 0; j < kernel[i].krnl.size(); j++) {
 419             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]);
 420         }
 421         std::fprintf(file, "\n");
 422     }
 423     std::fclose(file);
 424
 425
 426
 427     /*
 428     FILE *file;
 429     file = std::fopen("spiral_dist.txt", "w+");
 430     for (int i = 0; i < spiral_dist.size(); i++) {
 431         for (int j = 0; j < spiral_dist[i].size(); j++) {
 432             std::fprintf(file, "%3.2Lf\n", spiral_dist[i][j]);
 433         }
 434         std::fprintf(file, "\n");
 435     }
 436     std::fclose(file);*/
 437 }
 438
 439 void
 440 Spirals::writeOutput()
 441 {
 442
 443 }
 444
 445 /*! \brief
 446  * The main function for the analysis template.
 447  */
 448 int
 449 main(int argc, char *argv[])
 450 {
 451     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 452 }