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
  42 using namespace gmx;
  43
  44 struct crcl {
  45     long double x;
  46     long double y;
  47     long double z;
  48     long double r;
  49 };
  50
  51 crcl return_crcl(long double x1, long double y1, long double z1,
  52                  long double x2, long double y2, long double z2,
  53                  long double x3, long double y3, long double z3) {
  54     long double c, b, a, z, y, x, r;
  55     long double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12;
  56
  57     c = (z3 - z1 * x3 / x1 - (z2 - z1 * x2 / x1) / (y2 - y1 * x2 / x1)) / (1 - x3 / x1 - (1 - x2 / x1) / (y2 - y1 * x2 / x1));
  58     //std::cout << c << " ";
  59     //std::cout << y1 * (x2 / x1) << " ";
  60     b = (y2 - y1 * x2 / x1) / (1 - x2 / x1 - 1 / c * (z2 - z1 * x2 / x1));
  61     //std::cout << b << " ";
  62     a = (x1) / (1 - y1 / b - z1 / c);
  63     //std::cout << a << " ";
  64
  65     a1 = 1 / a;
  66     a2 = 1 / b;
  67     a3 = 1 / c;
  68     a4 = 1;
  69
  70     a5 = 2 * (x3 - x1);
  71     a6 = 2 * (y3 - y1);
  72     a7 = 2 * (z3 - z1);
  73     a8 = (x3 * x3 - x1 * x1) + (y3 * y3 - y1 * y1) + (z3 * z3 - z1 * z1);
  74
  75     a9  = 2 * (x3 - x2);
  76     a10 = 2 * (y3 - y2);
  77     a11 = 2 * (z3 - z2);
  78     a12 = (x3 * x3 - x2 * x2) + (y3 * y3 - y2 * y2) + (z3 * z3 - z2 * z2);
  79
  80     /*std::cout << a1 << " " << a2 << " " << a3 << " " << a4 << " "
  81                  << a5 << " " << a6 << " " << a7 << " " << a8 << " "
  82                     << a9 << " " << a10 << " " << a11 << " " << a12 << " ";*/
  83
  84     z = (a12 - a9 * a4 / a1 - (a10 - a9 * a2 / a1) * (a8 - a5 * a4 / a1) / (a6 - a5 * a2 / a1)) /
  85             (a11 - a9 * a3 / a1 - (a10 - a9 * a2 / a1) * (a7 - a5 * a3 / a1) / (a6 - a5 * a2 / a1));
  86     //std::cout << z << " ";
  87     y = 1 / (a6 - a5 * a2 / a1) * (a8 - a5 * a4 / a1 - z * (a7 - a5 * a3 / a1));
  88     //std::cout << y << " ";
  89     x = 1 / a1 * (a4 - a2 * y - a3 * z);
  90     //std::cout << x << " ";
  91     r = std::sqrt((x - x3) * (x - x3) + (y - y3) * (y - y3) + (z - z3) * (z - z3));
  92     //std::cout << r << " ";
  93
  94     crcl rtrn;
  95
  96     rtrn.x = x;
  97     rtrn.y = y;
  98     rtrn.z = z;
  99     rtrn.r = r;
 100
 101     return rtrn;
 102 }
 103
 104 /*! \brief
 105  * Template class to serve as a basis for user analysis tools.
 106  */
 107 class Spirals : public TrajectoryAnalysisModule
 108 {
 109     public:
 110         Spirals();
 111
 112         virtual void initOptions(IOptionsContainer          *options,
 113                                  TrajectoryAnalysisSettings *settings);
 114         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 115                                   const TopologyInformation        &top);
 116
 117         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 118                                   TrajectoryAnalysisModuleData *pdata);
 119
 120         virtual void finishAnalysis(int nframes);
 121         virtual void writeOutput();
 122
 123     private:
 124         class ModuleData;
 125
 126         std::string                      fnNdx_;
 127         double                           cutoff_;
 128         Selection                        refsel_;
 129
 130         AnalysisNeighborhood             nb_;
 131
 132         AnalysisData                     data_;
 133         AnalysisDataAverageModulePointer avem_;
 134
 135
 136         Selection                                           selec;
 137         std::vector< std::vector < RVec > >                 trajectory;
 138         std::vector< int >                                  index;
 139         int                                                 frames       = 0;
 140 };
 141
 142 Spirals::Spirals()
 143     : cutoff_(0.0)
 144 {
 145     registerAnalysisDataset(&data_, "avedist");
 146 }
 147
 148 void
 149 Spirals::initOptions(IOptionsContainer          *options,
 150                               TrajectoryAnalysisSettings *settings)
 151 {
 152     static const char *const desc[] = {
 153         "Analysis tool for finding molecular core."
 154     };
 155
 156     // Add the descriptive text (program help text) to the options
 157     settings->setHelpText(desc);
 158     // Add option for selecting a subset of atoms
 159     options->addOption(SelectionOption("select")
 160                            .store(&selec).required()
 161                            .description("Atoms that are considered as part of the excluded volume"));
 162     // Add option for output file name
 163     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 164                             .store(&fnNdx_).defaultBasename("rcore")
 165                             .description("Index file from the rcore"));
 166 }
 167
 168 // -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'
 169
 170 void
 171 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 172                                const TopologyInformation         & /*top*/)
 173 {
 174     index.resize(0);
 175     ConstArrayRef< int > atomind = selec.atomIndices();
 176     for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++)
 177     {
 178         index.push_back(*ai);
 179     }
 180 }
 181
 182 void
 183 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 184                                TrajectoryAnalysisModuleData *pdata)
 185 {
 186     trajectory.resize(frames + 1);
 187     trajectory[frames].resize(index.size());
 188     for (int i = 0; i < index.size(); i++)
 189     {
 190         trajectory[frames][i] = fr.x[index[i]];
 191     }
 192     frames++;
 193 }
 194
 195 void
 196 Spirals::finishAnalysis(int /*nframes*/)
 197 {
 198     long double x1, y1, z1, x2, y2, z2, x3, y3, z3;
 199     crcl test;
 200
 201     for (long double r = 1; r < 50; r++) {
 202         for (long double x0 = -5; x0 < 6; x0+=0.99) {
 203             for (long double y0 = -5; y0 < 6; y0+=0.99) {
 204                 x1 = x0 + r - r/100;
 205                 x2 = x0 + r/125;
 206                 x3 = x0 - r/7;
 207                 z1 = 1.01;
 208                 z2 = 1;
 209                 z3 = 0.99;
 210                 y1 = std::sqrt(r * r - (x1 - x0) * (x1 - x0)) + y0;
 211                 y2 = std::sqrt(r * r - (x2 - x0) * (x2 - x0)) + y0;
 212                 y3 = std::sqrt(r * r - (x3 - x0) * (x3 - x0)) + y0;
 213
 214                 /*std::cout << x1 << " " << y1 << " " << z1 << " "
 215                              << x2 << " " << y2 << " " << z2 << " "
 216                                 << x3 << " " << y3 << " " << z3 << "\n";*/
 217
 218                 test = return_crcl(x1, y1, z1, x2, y2, z2, x3, y3, z3);
 219
 220                 //std::cout.precision(1);
 221                 //std::cout.width(3);
 222                 //std::cout << std::abs(x0 - test.x) << " " << std::abs(y0 - test.y) << " " << std::abs(1 - test.z) << " " << std::abs(r - test.r) << "\n";
 223             }
 224         }
 225     }
 226
 227     for (int i = 0; i < 1; i++) {
 228         for (int j = 0; j < index.size() - 2; j++) {
 229             test = return_crcl( trajectory[i][j][0], trajectory[i][j][1], trajectory[i][j][2],
 230                                 trajectory[i][j + 1][0], trajectory[i][j + 1][1], trajectory[i][j + 1][2],
 231                                 trajectory[i][j + 2][0], trajectory[i][j + 2][1], trajectory[i][j + 2][2]);
 232             std::cout << test.x << " " << test.y << " " << test.z << " " << test.r << "\n";
 233         }
 234     }
 235
 236 }
 237
 238 void
 239 Spirals::writeOutput()
 240 {
 241
 242 }
 243
 244 /*! \brief
 245  * The main function for the analysis template.
 246  */
 247 int
 248 main(int argc, char *argv[])
 249 {
 250     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 251 }