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 #include <algorithm>
  39
  40 #include <gromacs/trajectoryanalysis.h>
  41 #include <gromacs/math/do_fit.h>
  42 #include <gromacs/utility/smalloc.h>
  43 #include "gromacs/selection/selection.h"
  44 #include "gromacs/selection/selectionoption.h"
  45
  46 using namespace gmx;
  47
  48 struct kernel_maxima {
  49     RVec x;
  50     RVec p;
  51     std::vector< RVec > krnl;
  52 };
  53
  54 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) {
  55     long double ret = 0;
  56     for (int i = 0; i < x.size(); i++) {
  57         ret +=
  58            sqrt (   pow (p2 * (x[i][2] - z0) - p3 * (x[i][1] - y0), 2) +
  59                         pow (p3 * (x[i][0] - x0) - p1 * (x[i][2] - z0), 2) +
  60                         pow (p1 * (x[i][1] - y0) - p2 * (x[i][0] - x0), 2)) /
  61            sqrt (p1 * p1 + p2 * p2 + p3 * p3);
  62     }
  63     return ret;
  64 }
  65
  66 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) {
  67     long double ret = 0;
  68     for (int i = 0; i < x.size(); i++) {
  69         ret +=
  70            (2 * p2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) + 2 * p3 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]))) /
  71            (2 *     sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  72                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  73                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
  74                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
  75     }
  76     return ret;
  77 }
  78
  79 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) {
  80     long double ret = 0;
  81     for (int i = 0; i < x.size(); i++) {
  82         ret +=
  83            -(2 * p1 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) - 2 * p3 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
  84             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  85                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  86                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
  87                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
  88     }
  89     return ret;
  90 }
  91
  92 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) {
  93     long double ret = 0;
  94     for (int i = 0; i < x.size(); i++) {
  95         ret +=
  96            -(2 * p1 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) + 2 * p2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
  97             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
  98                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
  99                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 100                     sqrt (p1 * p1 + p2 * p2 + p3 * p3));
 101     }
 102     return ret;
 103 }
 104
 105 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) {
 106     long double ret = 0;
 107     for (int i = 0; i < x.size(); i++) {
 108         ret +=
 109            -(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])) /
 110             (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 111                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 112                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 113                     sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 114             (p1 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 115                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 116                             pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 117                     pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 118     }
 119     return ret;
 120 }
 121
 122 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) {
 123     long double ret = 0;
 124     for (int i = 0; i < x.size(); i++) {
 125         ret +=
 126            (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])) /
 127            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 128                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 129                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 130                    sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 131            (p2 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 132                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 133                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 134                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 135     }
 136     return ret;
 137 }
 138
 139 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) {
 140     long double ret = 0;
 141     for (int i = 0; i < x.size(); i++) {
 142         ret +=
 143            (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])) /
 144            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 145                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 146                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
 147                    sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
 148            (p3 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
 149                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
 150                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
 151                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
 152     }
 153     return ret;
 154 }
 155
 156 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) {
 157     long double Ftemp = 0, FX = 0, FX0 = 0, FY0 = 0, FZ0 = 0, FP1 = 0, FP2 = 0, FP3 = 0, L0 = 0;
 158     while (true) {
 159         FX = Fx(x0, y0, z0, p1, p2, p3, x);
 160         FX0 = fx0(x0, y0, z0, p1, p2, p3, x);
 161         FY0 = fy0(x0, y0, z0, p1, p2, p3, x);
 162         FZ0 = fz0(x0, y0, z0, p1, p2, p3, x);
 163         FP1 = fp1(x0, y0, z0, p1, p2, p3, x);
 164         FP2 = fp2(x0, y0, z0, p1, p2, p3, x);
 165         FP3 = fp3(x0, y0, z0, p1, p2, p3, x);
 166
 167         L0 = 1;
 168         while (true) {
 169             Ftemp = Fx(x0 - L0 * FX0, y0 - L0 * FY0, z0 - L0 * FZ0, p1 - L0 * FP1, p2 - L0 * FP2, p3 - L0 * FP3, x);
 170             if (Ftemp - FX > 0) {
 171                 L0 /= 2;
 172             } else {
 173                 x0 -= L0 * FX0;
 174                 y0 -= L0 * FY0;
 175                 z0 -= L0 * FZ0;
 176                 p1 -= L0 * FP1;
 177                 p2 -= L0 * FP2;
 178                 p3 -= L0 * FP3;
 179                 if (FX - Ftemp < epsi) {
 180                     L0 = 0;
 181                 }
 182                 break;
 183             }
 184         }
 185
 186         if (L0 == 0) {
 187             break;
 188         }
 189     }
 190 }
 191
 192 RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p3, RVec x) {
 193     double lambda = (p1 * (x[0] - x0) + p2 * (x[1] - y0) + p3 * (x[2] - z0)) / (p1 * p1 + p2 * p2 + p3 * p3);
 194     RVec pro;
 195     pro[0] = x0 + p1 * lambda;
 196     pro[1] = y0 + p2 * lambda;
 197     pro[2] = z0 + p3 * lambda;
 198     return pro;
 199 }
 200
 201 void make_kernel (std::vector< RVec > temp, double epsi, std::vector< kernel_maxima > &kernel) {
 202     RVec mid, arrow;
 203
 204     mid[0] = 0;
 205     mid[1] = 0;
 206     mid[2] = 0;
 207
 208     arrow[0] = 0;
 209     arrow[1] = 0;
 210     arrow[2] = 0;
 211
 212     for (int i = 0; i < temp.size(); i++) {
 213         rvec_inc(temp[i], mid);
 214     }
 215     mid[0] /= temp.size();
 216     mid[1] /= temp.size();
 217     mid[2] /= temp.size();
 218     rvec_sub(temp.back(), temp.front(), arrow);
 219
 220     long double t1, t2, t3, t4, t5, t6;
 221     t1 = mid[0];
 222     t2 = mid[1];
 223     t3 = mid[2];
 224     t4 = arrow[0];
 225     t5 = arrow[1];
 226     t6 = arrow[2];
 227
 228     linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
 229
 230     mid[0] = t1;
 231     mid[1] = t2;
 232     mid[2] = t3;
 233     arrow[0] = t4;
 234     arrow[1] = t5;
 235     arrow[2] = t6;
 236
 237     kernel.back().x = mid;
 238     kernel.back().p = arrow;
 239     for (int i = 0; i < temp.size(); i++) {
 240         kernel.back().krnl.push_back(kernel_pro(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp[i]));
 241     }
 242 }
 243
 244 double left_right_turn (RVec a, RVec b, RVec c) {
 245     return  a[0] * b[1] * c[2] +
 246             a[1] * b[2] * c[0] +
 247             b[0] * c[1] * a[2] -
 248             a[2] * b[1] * c[0] -
 249             a[0] * b[2] * c[1] -
 250             a[1] * b[0] * c[2];
 251 }
 252
 253 void make_circles (std::vector< std::vector< std::vector< int > > > &circles, std::vector< RVec > temp, std::vector< kernel_maxima > kernel)  {
 254     bool st1 = true, st2 = false;
 255     double turn = -1, tempt;
 256     RVec a, b, c;
 257     rvec_sub(temp[0], kernel.back().krnl.front(), a);
 258     rvec_sub(kernel.back().krnl.front(), kernel.back().krnl.back(), b);
 259     for (int i = 1; i < temp.size(); i++) {
 260         rvec_sub(temp[i], kernel.back().krnl[i], c);
 261         tempt = left_right_turn(a, b, c);
 262         if (tempt > turn) {
 263             st1 = true;
 264         } else {
 265             st1 = false;
 266             st2 = true;
 267         }
 268         turn = tempt;
 269         if (st1 && !st2 || !st1 && st2) {
 270             if (circles.back().size() == 0) {
 271                 circles.back().resize(1);
 272             }
 273             circles.back().back().push_back(i);
 274         } else {
 275             circles.back().resize(circles.back().size() + 1);
 276             circles.back().back().push_back(i);
 277             st2 = false;
 278         }
 279     }
 280 }
 281
 282 /*! \brief
 283  * Template class to serve as a basis for user analysis tools.
 284  */
 285 class Spirals : public TrajectoryAnalysisModule
 286 {
 287     public:
 288         Spirals();
 289
 290         virtual void initOptions(IOptionsContainer          *options,
 291                                  TrajectoryAnalysisSettings *settings);
 292         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 293                                   const TopologyInformation        &top);
 294
 295         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 296                                   TrajectoryAnalysisModuleData *pdata);
 297
 298         virtual void finishAnalysis(int nframes);
 299         virtual void writeOutput();
 300
 301     private:
 302         class ModuleData;
 303
 304         std::string                      fnNdx_;
 305         double                           cutoff_;
 306         Selection                        refsel_;
 307
 308         AnalysisNeighborhood             nb_;
 309
 310         AnalysisData                     data_;
 311         AnalysisDataAverageModulePointer avem_;
 312
 313         SelectionList                           sel_;
 314         int                                     frames      = 0;
 315         double                                  epsi        = 0.00001;
 316         int                                     window      = 6;
 317
 318         std::vector< std::vector< RVec > >                  monomers;
 319         std::vector< kernel_maxima >                        kernel;
 320         std::vector< std::vector< std::vector< int > > >    circles;
 321         std::vector< std::vector< std::vector< int > > >    groups;
 322
 323         std::vector< std::vector< kernel_maxima > >                     window_kernel;
 324         std::vector< std::vector< std::vector< std::vector< int > > > > window_circles;
 325
 326 };
 327
 328 Spirals::Spirals()
 329     : cutoff_(0.0)
 330 {
 331     registerAnalysisDataset(&data_, "avedist");
 332 }
 333
 334 void
 335 Spirals::initOptions(IOptionsContainer          *options,
 336                               TrajectoryAnalysisSettings *settings)
 337 {
 338     static const char *const desc[] = {
 339         "Analysis tool for finding molecular core."
 340     };
 341
 342     // Add the descriptive text (program help text) to the options
 343     settings->setHelpText(desc);
 344     // Add option for output file name
 345     /*options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 346                             .store(&fnNdx_).defaultBasename("rcore")
 347                             .description("Index file from the rcore"));*/
 348     // Add option for window length constant
 349     options->addOption(gmx::IntegerOption("window_length")
 350                             .store(&window)
 351                             .description("window length for local parameters"));
 352     // Add option for selection list
 353     options->addOption(SelectionOption("select").storeVector(&sel_)
 354                            .required().dynamicMask().multiValue()
 355                            .description("Position to calculate distances for"));
 356 }
 357
 358 void
 359 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 360                                const TopologyInformation         & /*top*/)
 361 {
 362     kernel.resize(0);
 363     circles.resize(0);
 364     monomers.resize(0);
 365     groups.resize(0);
 366
 367
 368     window_kernel.resize(0);
 369     window_circles.resize(0);
 370 }
 371
 372 void
 373 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 374                                TrajectoryAnalysisModuleData *pdata)
 375 {
 376     const SelectionList &sel = pdata->parallelSelections(sel_);
 377     std::vector< RVec > temp;
 378         temp.resize(sel_.size());
 379     for (int i = 0; i < sel.size(); i++) {
 380         copy_rvec(sel[i].position(0).x(), temp[i]);
 381     }
 382
 383     int window_temp2 = sel_.size() - window + 1;
 384     if (window_kernel.size() == 0) {
 385         window_circles.resize(window_temp2);
 386         window_kernel.resize(window_temp2);
 387         for (int i = 0; i < window_temp2; i++) {
 388             window_circles[i].resize(0);
 389             window_kernel[i].resize(0);
 390         }
 391     }
 392
 393     monomers.resize(monomers.size() + 1);
 394     for (int i = 0; i < sel.size(); i++) {
 395         monomers.back().push_back(temp[i]);
 396     }
 397
 398     kernel.resize(kernel.size() + 1);
 399     make_kernel(temp, epsi, kernel);
 400
 401     circles.resize(circles.size() + 1);
 402     make_circles(circles, temp, kernel);
 403
 404     std::cout << "global parameters - end\n";
 405
 406     std::vector< std::vector< RVec > > window_temp;
 407     window_temp.resize(window_temp2);
 408     for (int i = 0; i < window_temp.size(); i++) {
 409         window_kernel[i].resize(window_kernel[i].size() + 1);
 410         window_kernel[i].back().p = kernel.back().p;
 411         window_kernel[i].back().x = kernel.back().x;
 412         window_temp[i].resize(window);
 413         for (int j = 0; j < window; j++) {
 414             window_temp[i].push_back(temp[i + j]);
 415             window_kernel[i].back().krnl.push_back(kernel.back().krnl[i + j]);
 416         }
 417     }
 418
 419     for (int i = 0; i < window_temp2; i++) {
 420         window_circles[i].resize(window_circles[i].size() + 1);
 421         make_circles(window_circles[i], window_temp[i], window_kernel[i]);
 422         std::cout << "local circles parameters " << i << " - end\n";
 423     }
 424
 425     /*groups.resize(groups.size() + 1);
 426     int itemp = circles.back().size() / 2;
 427     for (int j = 0; j < circles.back()[itemp].size(); j++) {
 428         groups.back().resize(groups.back().size() + 1);
 429         for (int i = itemp; i > -1; i--) {
 430             if (circles.back()[i].size() >= j) {
 431                 groups.back().back().push_back(circles.back()[i][circles.back()[i].size() - 1 - j]);
 432             }
 433         }
 434         for (int i = itemp + 1; i < circles.back().size(); i++) {
 435             if (circles.back()[i].size() >= circles.back()[itemp].size() - j) {
 436                 groups.back().back().push_back(circles.back()[i][circles.back()[i].size() - 1 - j]);
 437             }
 438         }
 439         std::sort(groups.back().begin(), groups.back().end());
 440     }*/
 441
 442     frames++;
 443 }
 444
 445 void
 446 Spirals::finishAnalysis(int /*nframes*/)
 447 {
 448     FILE *file;
 449
 450     file = std::fopen("linear_kernel.txt", "w+");
 451     for (int i = 0; i < kernel.size(); i++) {
 452         for (int j = 0; j < monomers[i].size(); j++) {
 453             std::fprintf(file, "%3.2f %3.2f %3.2f\n", monomers[i][j][0], monomers[i][j][1], monomers[i][j][2]);
 454         }
 455         std::fprintf(file, "\n");
 456         for (int j = 0; j < kernel[i].krnl.size(); j++) {
 457             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]);
 458         }
 459         std::fprintf(file, "\n\n");
 460     }
 461     std::fclose(file);
 462
 463     file = std::fopen("circles_points.txt", "w+");
 464     for (int i = 0; i < circles.size(); i++) {
 465         for (int j = 0; j < circles[i].size(); j++) {
 466             for (int k = 0; k < circles[i][j].size(); k++) {
 467                 std::fprintf(file, "%3.2d ", circles[i][j][k]);
 468             }
 469             std::fprintf(file, "\n");
 470         }
 471         std::fprintf(file, "\n");
 472     }
 473     std::fclose(file);
 474
 475     file = std::fopen("steps_Rspiral_Nmonomers.txt", "w+");
 476     for (int i = 0; i < circles.size(); i++) {
 477         std::fprintf(file, "Frame # %6.2d\n", i);
 478         std::fprintf(file, "Spiral steps:\n");
 479         for (int j = 0; j < circles[i].size(); j++) {
 480             std::fprintf(file, "%3.2f ", std::sqrt(distance2(kernel[i].krnl[circles[i][j].front() - 1], kernel[i].krnl[circles[i][j].back() - 1])));
 481         }
 482         std::fprintf(file, "\n");
 483         std::fprintf(file, "Spiral radii\n");
 484         for (int j = 0; j < circles[i].size(); j++) {
 485             float temp = 0;
 486             for (int k = 0; k < circles[i][j].size(); k++) {
 487                 temp += std::sqrt(distance2(kernel[i].krnl[circles[i][j][k] - 1], monomers[i][circles[i][j][k] - 1]));
 488                 std::fprintf(file, "%3.2f ", std::sqrt(distance2(kernel[i].krnl[circles[i][j][k] - 1], monomers[i][circles[i][j][k] - 1])));
 489             }
 490             std::fprintf(file, "average: %3.2f\n", temp / circles[i][j].size());
 491         }
 492         std::fprintf(file, "# of monomers per coil:\n");
 493         for (int j = 0; j < circles[i].size(); j++) {
 494             std::fprintf(file, "%3.2d ", circles[i][j].size());
 495         }
 496         std::fprintf(file, "\n\n");
 497     }
 498     std::fclose(file);
 499
 500     file = std::fopen("LocalSteps_Rspiral_Nmonomers.txt", "w+");
 501     for (int m = 0; m < window_circles.front().size(); m++) {
 502         for (int i = 0; i < window_circles.size(); i++) {
 503             std::fprintf(file, "Frame # %6.2d | window # %3.2d\n", m, i);
 504             std::fprintf(file, "Spiral steps:\n");
 505             for (int j = 0; j < window_circles[i][m].size(); j++) {
 506                 std::fprintf(file, "%3.2f ", std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j].front() - 1], window_kernel[i][m].krnl[window_circles[i][m][j].back() - 1])));
 507             }
 508             std::cout << "spiral steps - ok\n";
 509             std::fprintf(file, "\n");
 510             std::fprintf(file, "Spiral radii\n");
 511             for (int j = 0; j < window_circles[i][m].size(); j++) {
 512                 float temp = 0;
 513                 for (int k = 0; k < window_circles[i][m][j].size(); k++) {
 514                     temp += std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j][k] - 1], monomers[i][window_circles[i][m][j][k] - 1 + i]));
 515                     std::fprintf(file, "%3.2f ", std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j][k] - 1], monomers[i][window_circles[i][m][j][k] - 1 + i])));
 516                 }
 517                 std::fprintf(file, "average: %3.2f\n", temp / window_circles[i][m][j].size());
 518             }
 519             std::cout << "spiral radii - ok\n";
 520             std::fprintf(file, "# of monomers per coil:\n");
 521             for (int j = 0; j < window_circles[i][m].size(); j++) {
 522                 std::fprintf(file, "%3.2d ", window_circles[i][m][j].size());
 523             }
 524             std::fprintf(file, "\n\n");
 525             std::cout << "spiral numbers - ok\n";
 526         }
 527     }
 528     std::fclose(file);
 529 }
 530
 531 void
 532 Spirals::writeOutput()
 533 {
 534
 535 }
 536
 537 /*! \brief
 538  * The main function for the analysis template.
 539  */
 540 int
 541 main(int argc, char *argv[])
 542 {
 543     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 544 }