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     circles.back().resize(1);
 260     circles.back().back().push_back(0);
 261     for (int i = 1; i < temp.size(); i++) {
 262         rvec_sub(temp[i], kernel.back().krnl[i], c);
 263         tempt = left_right_turn(a, b, c);
 264         if (tempt > turn) {
 265             st1 = true;
 266         } else {
 267             st1 = false;
 268             st2 = true;
 269         }
 270         turn = tempt;
 271         if (st1 && !st2 || !st1 && st2) {
 272             if (circles.back().size() == 0) {
 273                 circles.back().resize(1);
 274             }
 275             circles.back().back().push_back(i);
 276         } else {
 277             circles.back().resize(circles.back().size() + 1);
 278             circles.back().back().push_back(i);
 279             st2 = false;
 280         }
 281     }
 282 }
 283
 284 /*! \brief
 285  * Template class to serve as a basis for user analysis tools.
 286  */
 287 class Spirals : public TrajectoryAnalysisModule
 288 {
 289     public:
 290         Spirals();
 291
 292         virtual void initOptions(IOptionsContainer          *options,
 293                                  TrajectoryAnalysisSettings *settings);
 294         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 295                                   const TopologyInformation        &top);
 296
 297         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 298                                   TrajectoryAnalysisModuleData *pdata);
 299
 300         virtual void finishAnalysis(int nframes);
 301         virtual void writeOutput();
 302
 303     private:
 304         class ModuleData;
 305
 306         std::string                      fnNdx_;
 307         double                           cutoff_;
 308         Selection                        refsel_;
 309
 310         AnalysisNeighborhood             nb_;
 311
 312         AnalysisData                     data_;
 313         AnalysisDataAverageModulePointer avem_;
 314
 315         SelectionList                           sel_;
 316         int                                     frames      = 0;
 317         double                                  epsi        = 0.00001;
 318         int                                     window      = 6;
 319
 320         std::vector< std::vector< RVec > >                  monomers;
 321         std::vector< kernel_maxima >                        kernel;
 322         std::vector< std::vector< std::vector< int > > >    circles;
 323         std::vector< std::vector< std::vector< int > > >    groups;
 324
 325         std::vector< std::vector< kernel_maxima > >                     window_kernel;
 326         std::vector< std::vector< std::vector< std::vector< int > > > > window_circles;
 327
 328 };
 329
 330 Spirals::Spirals()
 331     : cutoff_(0.0)
 332 {
 333     registerAnalysisDataset(&data_, "avedist");
 334 }
 335
 336 void
 337 Spirals::initOptions(IOptionsContainer          *options,
 338                               TrajectoryAnalysisSettings *settings)
 339 {
 340     static const char *const desc[] = {
 341         "Analysis tool for finding molecular core."
 342     };
 343
 344     // Add the descriptive text (program help text) to the options
 345     settings->setHelpText(desc);
 346     // Add option for output file name
 347     /*options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 348                             .store(&fnNdx_).defaultBasename("rcore")
 349                             .description("Index file from the rcore"));*/
 350     // Add option for window length constant
 351     options->addOption(gmx::IntegerOption("window_length")
 352                             .store(&window)
 353                             .description("window length for local parameters"));
 354     // Add option for selection list
 355     options->addOption(SelectionOption("select").storeVector(&sel_)
 356                            .required().dynamicMask().multiValue()
 357                            .description("Position to calculate distances for"));
 358 }
 359
 360 void
 361 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
 362                                const TopologyInformation         & /*top*/)
 363 {
 364     kernel.resize(0);
 365     circles.resize(0);
 366     monomers.resize(0);
 367     groups.resize(0);
 368
 369
 370     window_kernel.resize(0);
 371     window_circles.resize(0);
 372 }
 373
 374 void
 375 Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 376                                TrajectoryAnalysisModuleData *pdata)
 377 {
 378     const SelectionList &sel = pdata->parallelSelections(sel_);
 379     std::vector< RVec > temp;
 380         temp.resize(sel_.size());
 381     for (int i = 0; i < sel.size(); i++) {
 382         copy_rvec(sel[i].position(0).x(), temp[i]);
 383     }
 384
 385     int window_temp2 = sel_.size() - window + 1;
 386     if (window_kernel.size() == 0) {
 387         window_circles.resize(window_temp2);
 388         window_kernel.resize(window_temp2);
 389         for (int i = 0; i < window_temp2; i++) {
 390             window_circles[i].resize(0);
 391             window_kernel[i].resize(0);
 392         }
 393     }
 394
 395     monomers.resize(monomers.size() + 1);
 396     for (int i = 0; i < sel.size(); i++) {
 397         monomers.back().push_back(temp[i]);
 398     }
 399
 400     kernel.resize(kernel.size() + 1);
 401     make_kernel(temp, epsi, kernel);
 402
 403     circles.resize(circles.size() + 1);
 404     make_circles(circles, temp, kernel);
 405
 406     std::cout << "global parameters - end\n";
 407
 408     std::vector< std::vector< RVec > > window_temp;
 409     window_temp.resize(window_temp2);
 410     for (int i = 0; i < window_temp2; i++) {
 411         window_kernel[i].resize(window_kernel[i].size() + 1);
 412         window_kernel[i].back().p = kernel.back().p;
 413         window_kernel[i].back().x = kernel.back().x;
 414         window_temp[i].resize(window);
 415         for (int j = 0; j < window; j++) {
 416             window_temp[i].push_back(temp[i + j]);
 417             window_kernel[i].back().krnl.push_back(kernel.back().krnl[i + j]);
 418         }
 419     }
 420
 421     for (int i = 0; i < window_temp2; i++) {
 422         window_circles[i].resize(window_circles[i].size() + 1);
 423         make_circles(window_circles[i], window_temp[i], window_kernel[i]);
 424         std::cout << window_temp[i].size() << " " << window_kernel[i].back().krnl.size() << "\n";
 425         for (int j = 0; j < window_circles[i].back().size(); j++) {
 426             std::cout << window_circles[i].back()[j].size() << " ";
 427         }
 428         std::cout << "\n";
 429     }
 430
 431     /*groups.resize(groups.size() + 1);
 432     int itemp = circles.back().size() / 2;
 433     for (int j = 0; j < circles.back()[itemp].size(); j++) {
 434         groups.back().resize(groups.back().size() + 1);
 435         for (int i = itemp; i > -1; i--) {
 436             if (circles.back()[i].size() >= j) {
 437                 groups.back().back().push_back(circles.back()[i][circles.back()[i].size() - 1 - j]);
 438             }
 439         }
 440         for (int i = itemp + 1; i < circles.back().size(); i++) {
 441             if (circles.back()[i].size() >= circles.back()[itemp].size() - j) {
 442                 groups.back().back().push_back(circles.back()[i][circles.back()[i].size() - 1 - j]);
 443             }
 444         }
 445         std::sort(groups.back().begin(), groups.back().end());
 446     }*/
 447
 448     frames++;
 449 }
 450
 451 void
 452 Spirals::finishAnalysis(int /*nframes*/)
 453 {
 454     FILE *file;
 455
 456     file = std::fopen("linear_kernel.txt", "w+");
 457     for (int i = 0; i < kernel.size(); i++) {
 458         for (int j = 0; j < monomers[i].size(); j++) {
 459             std::fprintf(file, "%3.2f %3.2f %3.2f\n", monomers[i][j][0], monomers[i][j][1], monomers[i][j][2]);
 460         }
 461         std::fprintf(file, "\n");
 462         for (int j = 0; j < kernel[i].krnl.size(); j++) {
 463             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]);
 464         }
 465         std::fprintf(file, "\n\n");
 466     }
 467     std::fclose(file);
 468
 469     file = std::fopen("circles_points.txt", "w+");
 470     for (int i = 0; i < circles.size(); i++) {
 471         for (int j = 0; j < circles[i].size(); j++) {
 472             for (int k = 0; k < circles[i][j].size(); k++) {
 473                 std::fprintf(file, "%3.2d ", circles[i][j][k]);
 474             }
 475             std::fprintf(file, "\n");
 476         }
 477         std::fprintf(file, "\n");
 478     }
 479     std::fclose(file);
 480
 481     file = std::fopen("steps_Rspiral_Nmonomers.txt", "w+");
 482     for (int i = 0; i < circles.size(); i++) {
 483         std::fprintf(file, "Frame # %6.2d\n", i);
 484         std::fprintf(file, "Spiral steps:\n");
 485         for (int j = 0; j < circles[i].size(); j++) {
 486             std::fprintf(file, "%3.2f ", std::sqrt(distance2(kernel[i].krnl[circles[i][j].front()], kernel[i].krnl[circles[i][j].back()])));
 487         }
 488         std::fprintf(file, "\n");
 489         std::fprintf(file, "Spiral radii\n");
 490         for (int j = 0; j < circles[i].size(); j++) {
 491             float temp = 0;
 492             for (int k = 0; k < circles[i][j].size(); k++) {
 493                 temp += std::sqrt(distance2(kernel[i].krnl[circles[i][j][k]], monomers[i][circles[i][j][k]]));
 494                 std::fprintf(file, "%3.2f ", std::sqrt(distance2(kernel[i].krnl[circles[i][j][k]], monomers[i][circles[i][j][k]])));
 495             }
 496             std::fprintf(file, "average: %3.2f\n", temp / circles[i][j].size());
 497         }
 498         std::fprintf(file, "# of monomers per coil:\n");
 499         for (int j = 0; j < circles[i].size(); j++) {
 500             std::fprintf(file, "%3.2lu ", circles[i][j].size());
 501         }
 502         std::fprintf(file, "\n\n");
 503     }
 504     std::fclose(file);
 505
 506     file = std::fopen("LocalSteps_Rspiral_Nmonomers.txt", "w+");
 507     for (int i = 0; i < window_circles.front().size(); i++) {
 508         for (int j = 0; j < window_circles.size(); j++) {
 509             for (int k = 0; k < window_circles[j][i].size(); k++) {
 510                 std::fprintf(file, "%3.2lu ", window_circles[j][i][k].size());
 511             }
 512             std::fprintf(file, "\n");
 513         }
 514     }
 515     std::fclose(file);
 516
 517     /*file = std::fopen("LocalSteps_Rspiral_Nmonomers.txt", "w+");
 518     std::cout << " " << window_circles.front().size() << "\n";
 519     for (int m = 0; m < window_circles.front().size(); m++) {
 520         for (int i = 0; i < window_circles.size(); i++) {
 521             std::fprintf(file, "Frame # %6.2d | window # %3.2d\n", m, i);
 522             std::fprintf(file, "Spiral steps:\n");
 523             for (int j = 0; j < window_circles[i][m].size(); j++) {
 524                 std::fprintf(file, "%3.2f ", std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j].front()], window_kernel[i][m].krnl[window_circles[i][m][j].back()])));
 525             }
 526             std::fprintf(file, "\n");
 527             std::fprintf(file, "Spiral radii\n");
 528             for (int j = 0; j < window_circles[i][m].size(); j++) {
 529                 float temp = 0;
 530                 for (int k = 0; k < window_circles[i][m][j].size(); k++) {
 531                     temp += std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j][k]], monomers[m][window_circles[i][m][j][k] + i]));
 532                     std::fprintf(file, "%3.2f ", std::sqrt(distance2(window_kernel[i][m].krnl[window_circles[i][m][j][k]], monomers[m][window_circles[i][m][j][k] + i])));
 533                 }
 534                 std::fprintf(file, "average: %3.2f\n", temp / window_circles[i][m][j].size());
 535             }
 536             std::fprintf(file, "# of monomers per coil:\n");
 537             for (int j = 0; j < window_circles[i][m].size(); j++) {
 538                 std::fprintf(file, "%3.2lu ", window_circles[i][m][j].size());
 539             }
 540             std::fprintf(file, "\n\n");
 541         }
 542     }
 543     std::fclose(file);*/
 544 }
 545
 546 void
 547 Spirals::writeOutput()
 548 {
 549
 550 }
 551
 552 /*! \brief
 553  * The main function for the analysis template.
 554  */
 555 int
 556 main(int argc, char *argv[])
 557 {
 558     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
 559 }