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