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
  38 #include <gromacs/trajectoryanalysis.h>
  39 #include <gromacs/math/do_fit.h>
  40 #include <gromacs/utility/smalloc.h>
  41 #include "gromacs/selection/selection.h"
  42 #include "gromacs/selection/selectionoption.h"
  43
  44 using namespace gmx;
  45
  46 struct crcl {
  47     long double x;
  48     long double y;
  49     long double z;
  50     long double r;
  51 };
  52
  53 crcl return_crcl(long double x1, long double y1, long double z1,
  54                  long double x2, long double y2, long double z2,
  55                  long double x3, long double y3, long double z3) {
  56     long double c, b, a, z, y, x, r;
  57     long double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12;
  58
  59     c = (z3 - z1 * x3 / x1 - (z2 - z1 * x2 / x1) / (y2 - y1 * x2 / x1)) / (1 - x3 / x1 - (1 - x2 / x1) / (y2 - y1 * x2 / x1));
  60     b = (y2 - y1 * x2 / x1) / (1 - x2 / x1 - 1 / c * (z2 - z1 * x2 / x1));
  61     a = (x1) / (1 - y1 / b - z1 / c);
  62
  63     a1 = 1 / a;
  64     a2 = 1 / b;
  65     a3 = 1 / c;
  66     a4 = 1;
  67
  68     a5 = 2 * (x3 - x1);
  69     a6 = 2 * (y3 - y1);
  70     a7 = 2 * (z3 - z1);
  71     a8 = (x3 * x3 - x1 * x1) + (y3 * y3 - y1 * y1) + (z3 * z3 - z1 * z1);
  72
  73     a9  = 2 * (x3 - x2);
  74     a10 = 2 * (y3 - y2);
  75     a11 = 2 * (z3 - z2);
  76     a12 = (x3 * x3 - x2 * x2) + (y3 * y3 - y2 * y2) + (z3 * z3 - z2 * z2);
  77
  78     z = (a12 - a9 * a4 / a1 - (a10 - a9 * a2 / a1) * (a8 - a5 * a4 / a1) / (a6 - a5 * a2 / a1)) /
  79             (a11 - a9 * a3 / a1 - (a10 - a9 * a2 / a1) * (a7 - a5 * a3 / a1) / (a6 - a5 * a2 / a1));
  80     y = 1 / (a6 - a5 * a2 / a1) * (a8 - a5 * a4 / a1 - z * (a7 - a5 * a3 / a1));
  81     x = 1 / a1 * (a4 - a2 * y - a3 * z);
  82     r = std::sqrt((x - x3) * (x - x3) + (y - y3) * (y - y3) + (z - z3) * (z - z3));
  83
  84     crcl rtrn;
  85
  86     rtrn.x = x;
  87     rtrn.y = y;
  88     rtrn.z = z;
  89     rtrn.r = r;
  90
  91     return rtrn;
  92 }
  93
  94 long double lr_turn(std::vector< long double > a, std::vector< long double > b, std::vector< long double > c) {
  95     return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + b[0] * c[1] * a[2] -
  96             (a[2] * b[1] * c[0] + a[0] * b[2] * c[1] + a[1] * b[0] * c[2]);
  97 }
  98
  99 long double kernel_dist(crcl a, crcl b) {
 100     return std::sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
 101 }
 102
 103 /*! \brief
 104  * Template class to serve as a basis for user analysis tools.
 105  */
 106 class Spirals : public TrajectoryAnalysisModule
 107 {
 108     public:
 109         Spirals();
 110
 111         virtual void initOptions(IOptionsContainer          *options,
 112                                  TrajectoryAnalysisSettings *settings);
 113         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 114                                   const TopologyInformation        &top);
 115
 116         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 117                                   TrajectoryAnalysisModuleData *pdata);
 118
 119         virtual void finishAnalysis(int nframes);
 120         virtual void writeOutput();
 121
 122     private:
 123         class ModuleData;
 124
 125         std::string                      fnNdx_;
 126         double                           cutoff_;
 127         Selection                        refsel_;
 128
 129         AnalysisNeighborhood             nb_;
 130
 131         AnalysisData                     data_;
 132         AnalysisDataAverageModulePointer avem_;
 133
 134         SelectionList                                       sel_;
 135         int                                                 frames       = 0;
 136         std::vector< std::vector< crcl > >                  kernel;
 137         std::vector< std::vector< std::vector< int > > >    circles;
 138         std::vector< std::vector< long double > >           spiral_dist;
 139 };
 140
 141 Spirals::Spirals()
 142     : cutoff_(0.0)
 143 {
 144     registerAnalysisDataset(&data_, "avedist");
 145 }
 146
 147 void
 148 Spirals::initOptions(IOptionsContainer          *options,
 149                               TrajectoryAnalysisSettings *settings)
 150 {
 151     static const char *const desc[] = {
 152         "Analysis tool for finding molecular core."
 153     };
 154
 155     // Add the descriptive text (program help text) to the options
 156     settings->setHelpText(desc);
 157     // Add option for output file name
 158     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 159                             .store(&fnNdx_).defaultBasename("rcore")
 160                             .description("Index file from the rcore"));
 161     // Add option for selection list
 162     options->addOption(SelectionOption("select").storeVector(&sel_)
 163                            .required().dynamicMask().multiValue()
 164                            .description("Position to calculate distances for"));
 165 }
 166
 167 // -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'
 168
 169 void
 170 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 171                                const TopologyInformation         & /*top*/)
 172 {
 173
 174 }
 175
 176 void
 177 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 178                                TrajectoryAnalysisModuleData *pdata)
 179 {
 180     kernel.resize(frames + 1);
 181     circles.resize(frames + 1);
 182     spiral_dist.resize(frames + 1);
 183
 184     const SelectionList &sel = pdata->parallelSelections(sel_);
 185     std::vector< RVec > temp;
 186     temp.resize(sel_.size());
 187     for (int i = 0; i < sel.size(); i++) {
 188         copy_rvec(sel[i].position(0).x(), temp[i]);
 189     }
 190
 191     // центры окружностей и их радиусы
 192     for (int i = 0; i < temp.size() - 2; i++) {
 193         kernel[frames].push_back(return_crcl(   temp[i][0],     temp[i][1],     temp[i][2],
 194                                         temp[i + 1][0], temp[i + 1][1], temp[i + 1][2],
 195                                         temp[i + 2][0], temp[i + 2][1], temp[i + 2][2]));
 196     }
 197
 198     // распределение точек по виткам
 199     std::vector< long double > a, b, c;
 200     a.resize(3);
 201     b.resize(3);
 202     c.resize(3);
 203
 204     a[0] = temp[0][0] - kernel[frames].begin()->x;
 205     a[1] = temp[0][1] - kernel[frames].begin()->y;
 206     a[2] = temp[0][2] - kernel[frames].begin()->z;
 207
 208     c[0] = kernel[frames].end()->x - kernel[frames].begin()->x;
 209     c[1] = kernel[frames].end()->y - kernel[frames].begin()->y;
 210     c[2] = kernel[frames].end()->z - kernel[frames].begin()->z;
 211
 212     std::vector< int > empty;
 213     empty.resize(0);
 214
 215
 216     circles[frames].push_back(empty);
 217     circles[frames][0].push_back(0);
 218
 219     long double sign = 0, prev_sign = 0;
 220     for (int i = 1; i < temp.size(); i++) {
 221         if (i < kernel[frames].size()) {
 222             b[0] = temp[i][0] - kernel[frames][i].x;
 223             b[1] = temp[i][1] - kernel[frames][i].y;
 224             b[2] = temp[i][2] - kernel[frames][i].z;
 225         } else {
 226             b[0] = temp[i][0] - kernel[frames].end()->x;
 227             b[1] = temp[i][1] - kernel[frames].end()->y;
 228             b[2] = temp[i][2] - kernel[frames].end()->z;
 229         }
 230         long double local_sign = lr_turn(a, b ,c);
 231         if (sign == 0) {
 232             sign = local_sign;
 233             prev_sign = local_sign;
 234             circles[frames].back().push_back(i);
 235         } else {
 236             if (std::signbit(local_sign) == std::signbit(prev_sign)) {
 237                 circles[frames].back().push_back(i);
 238             } else {
 239                 if (std::signbit(local_sign) != std::signbit(sign)) {
 240                     circles[frames].back().push_back(i);
 241                 } else {
 242                     circles[frames].push_back(empty);
 243                     circles[frames].back().push_back(i);
 244                 }
 245             }
 246             prev_sign = local_sign;
 247         }
 248     }
 249
 250     // шаг спирали
 251     spiral_dist[frames].resize(circles[frames].size(), 0);
 252     for (int i = 0; i < circles[frames].size(); i++) {
 253         if (i == 0) {
 254             spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]) +
 255                              kernel_dist(kernel[frames][circles[frames][i].back()], kernel[frames][circles[frames][i + 1].front()]) / 2;
 256         } else if (i != circles[frames].size() - 1) {
 257             spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i - 1].back()], kernel[frames][circles[frames][i].front()]) / 2 +
 258                              kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]) +
 259                              kernel_dist(kernel[frames][circles[frames][i].back()], kernel[frames][circles[frames][i + 1].front()]) / 2;
 260         } else if (i == circles[frames].size() - 1) {
 261             spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i - 1].back()], kernel[frames][circles[frames][i].front()]) / 2 +
 262                              kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]);
 263         }
 264     }
 265
 266    frames++;
 267 }
 268
 269 void
 270 Spirals::finishAnalysis(int /*nframes*/)
 271 {
 272
 273 }
 274
 275 void
 276 Spirals::writeOutput()
 277 {
 278
 279 }
 280
 281 /*! \brief
 282  * The main function for the analysis template.
 283  */
 284 int
 285 main(int argc, char *argv[])
 286 {
 287     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 288 }