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