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         SelectionList                                       sel_;
 138         std::vector< std::vector < RVec > >                 trajectory;
 139         std::vector< std::vector< int > >                   index_all;
 140         std::vector< int >                                  index;
 141         int                                                 frames       = 0;
 142 };
 143
 144 Spirals::Spirals()
 145     : cutoff_(0.0)
 146 {
 147     registerAnalysisDataset(&data_, "avedist");
 148 }
 149
 150 void
 151 Spirals::initOptions(IOptionsContainer          *options,
 152                               TrajectoryAnalysisSettings *settings)
 153 {
 154     static const char *const desc[] = {
 155         "Analysis tool for finding molecular core."
 156     };
 157
 158     // Add the descriptive text (program help text) to the options
 159     settings->setHelpText(desc);
 160     // Add option for output file name
 161     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 162                             .store(&fnNdx_).defaultBasename("rcore")
 163                             .description("Index file from the rcore"));
 164     // Add option for selection list
 165     options->addOption(SelectionOption("Select groups").storeVector(&sel_)
 166                            .required().dynamicMask().multiValue()
 167                            .description("Domains to form rigid skeleton"));
 168 }
 169
 170 // -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'
 171
 172 void
 173 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 174                                const TopologyInformation         & /*top*/)
 175 {
 176     std::vector< int > a;
 177     a.resize(0);
 178     index_all.resize(0);
 179     index.resize(0);
 180     for (int i = 0; i < sel_.size(); i++) {
 181         index_all.push_back(a);
 182         ConstArrayRef< int > atomind = sel_[i].atomIndices();
 183         for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 184             index_all.back().push_back(*ai);
 185             index.push_back(*ai);
 186         }
 187     }
 188 }
 189
 190 void
 191 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 192                                TrajectoryAnalysisModuleData *pdata)
 193 {
 194     trajectory.resize(frames + 1);
 195     trajectory[frames].resize(index.size());
 196     for (int i = 0; i < index.size(); i++)
 197     {
 198         trajectory[frames][i] = fr.x[index[i]];
 199     }
 200     frames++;
 201 }
 202
 203 void
 204 Spirals::finishAnalysis(int /*nframes*/)
 205 {
 206     long double x1, y1, z1, x2, y2, z2, x3, y3, z3;
 207     crcl test;
 208
 209     for (long double r = 1; r < 50; r++) {
 210         for (long double x0 = -5; x0 < 6; x0+=0.99) {
 211             for (long double y0 = -5; y0 < 6; y0+=0.99) {
 212                 x1 = x0 + r - r/100;
 213                 x2 = x0 + r/125;
 214                 x3 = x0 - r/7;
 215                 z1 = 1.01;
 216                 z2 = 1;
 217                 z3 = 0.99;
 218                 y1 = std::sqrt(r * r - (x1 - x0) * (x1 - x0)) + y0;
 219                 y2 = std::sqrt(r * r - (x2 - x0) * (x2 - x0)) + y0;
 220                 y3 = std::sqrt(r * r - (x3 - x0) * (x3 - x0)) + y0;
 221
 222                 /*std::cout << x1 << " " << y1 << " " << z1 << " "
 223                              << x2 << " " << y2 << " " << z2 << " "
 224                                 << x3 << " " << y3 << " " << z3 << "\n";*/
 225
 226                 test = return_crcl(x1, y1, z1, x2, y2, z2, x3, y3, z3);
 227
 228                 //std::cout.precision(1);
 229                 //std::cout.width(3);
 230                 //std::cout << std::abs(x0 - test.x) << " " << std::abs(y0 - test.y) << " " << std::abs(1 - test.z) << " " << std::abs(r - test.r) << "\n";
 231             }
 232         }
 233     }
 234
 235     std::vector< std::vector< long double > > masses;
 236     masses.resize(0);
 237     int count = 0;
 238     for (int i = 0; i < index_all.size(); i++) {
 239         std::vector< long double > center;
 240         center.resize(0, 3);
 241         for (int j = 0; j < index_all[i].size(); j++) {
 242             center[0] += trajectory[0][count][0];
 243             center[1] += trajectory[0][count][1];
 244             center[2] += trajectory[0][count][2];
 245             count++;
 246         }
 247         center[0] /= index_all[i].size();
 248         center[1] /= index_all[i].size();
 249         center[2] /= index_all[i].size();
 250         masses.push_back(center);
 251     }
 252
 253
 254     for (int i = 0; i < 1; i++) {
 255         for (int j = 0; j < masses.size() - 2; j++) {
 256             /*test = return_crcl( trajectory[i][j][0], trajectory[i][j][1], trajectory[i][j][2],
 257                                 trajectory[i][j + 1][0], trajectory[i][j + 1][1], trajectory[i][j + 1][2],
 258                                 trajectory[i][j + 2][0], trajectory[i][j + 2][1], trajectory[i][j + 2][2]);*/
 259             test = return_crcl( masses[j][0], masses[j][1], masses[j][2],
 260                                 masses[j + 1][0], masses[j + 1][1], masses[j + 1][2],
 261                                 masses[j + 2][0], masses[j + 2][1], masses[j + 2][2]);
 262             std::cout << test.x << " " << test.y << " " << test.z << " " << test.r << "\n";
 263         }
 264     }
 265
 266 }
 267
 268 void
 269 Spirals::writeOutput()
 270 {
 271
 272 }
 273
 274 /*! \brief
 275  * The main function for the analysis template.
 276  */
 277 int
 278 main(int argc, char *argv[])
 279 {
 280     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 281 }