src/gromacs/gmxana/gmx_order.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/cmat.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/gmxana/gstat.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/pbcutil/rmpbc.h"
  58 #include "gromacs/topology/index.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/trajectory/trajectoryframe.h"
  61 #include "gromacs/utility/arraysize.h"
  62 #include "gromacs/utility/cstringutil.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/gmxassert.h"
  66 #include "gromacs/utility/smalloc.h"
  67
  68 /****************************************************************************/
  69 /* This program calculates the order parameter per atom for an interface or */
  70 /* bilayer, averaged over time.                                             */
  71 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij))                                 */
  72 /* It is assumed that the order parameter with respect to a box-axis        */
  73 /* is calculated. In that case i = j = axis, and delta(ij) = 1.             */
  74 /*                                                                          */
  75 /* Peter Tieleman,  April 1995                                              */
  76 /* P.J. van Maaren, November 2005     Added tetrahedral stuff               */
  77 /****************************************************************************/
  78
  79 static void find_nearest_neighbours(PbcType     pbcType,
  80                                     int         natoms,
  81                                     matrix      box,
  82                                     rvec        x[],
  83                                     int         maxidx,
  84                                     const int   index[],
  85                                     real*       sgmean,
  86                                     real*       skmean,
  87                                     int         nslice,
  88                                     int         slice_dim,
  89                                     real        sgslice[],
  90                                     real        skslice[],
  91                                     gmx_rmpbc_t gpbc)
  92 {
  93     int   ix, jx, nsgbin, *sgbin;
  94     int   i, ibin, j, k, l, n, *nn[4];
  95     rvec  dx, rj, rk, urk, urj;
  96     real  cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
  97     t_pbc pbc;
  98     int   sl_index;
  99     real* sl_count;
 100     real  onethird = 1.0 / 3.0;
 101     /*  dmat = init_mat(maxidx, FALSE); */
 102     box2 = box[XX][XX] * box[XX][XX];
 103     snew(sl_count, nslice);
 104     for (i = 0; (i < 4); i++)
 105     {
 106         snew(r_nn[i], natoms);
 107         snew(nn[i], natoms);
 108
 109         for (j = 0; (j < natoms); j++)
 110         {
 111             r_nn[i][j] = box2;
 112         }
 113     }
 114
 115     snew(sgmol, maxidx);
 116     snew(skmol, maxidx);
 117
 118     /* Must init pbc every step because of pressure coupling */
 119     set_pbc(&pbc, pbcType, box);
 120
 121     gmx_rmpbc(gpbc, natoms, box, x);
 122
 123     nsgbin = 2001; // Calculated as (1 + 1/0.0005)
 124     snew(sgbin, nsgbin);
 125
 126     *sgmean = 0.0;
 127     *skmean = 0.0;
 128     l       = 0;
 129     for (i = 0; (i < maxidx); i++) /* loop over index file */
 130     {
 131         ix = index[i];
 132         for (j = 0; (j < maxidx); j++)
 133         {
 134             if (i == j)
 135             {
 136                 continue;
 137             }
 138
 139             jx = index[j];
 140
 141             pbc_dx(&pbc, x[ix], x[jx], dx);
 142             r2 = iprod(dx, dx);
 143
 144             /* set_mat_entry(dmat,i,j,r2); */
 145
 146             /* determine the nearest neighbours */
 147             if (r2 < r_nn[0][i])
 148             {
 149                 r_nn[3][i] = r_nn[2][i];
 150                 nn[3][i]   = nn[2][i];
 151                 r_nn[2][i] = r_nn[1][i];
 152                 nn[2][i]   = nn[1][i];
 153                 r_nn[1][i] = r_nn[0][i];
 154                 nn[1][i]   = nn[0][i];
 155                 r_nn[0][i] = r2;
 156                 nn[0][i]   = j;
 157             }
 158             else if (r2 < r_nn[1][i])
 159             {
 160                 r_nn[3][i] = r_nn[2][i];
 161                 nn[3][i]   = nn[2][i];
 162                 r_nn[2][i] = r_nn[1][i];
 163                 nn[2][i]   = nn[1][i];
 164                 r_nn[1][i] = r2;
 165                 nn[1][i]   = j;
 166             }
 167             else if (r2 < r_nn[2][i])
 168             {
 169                 r_nn[3][i] = r_nn[2][i];
 170                 nn[3][i]   = nn[2][i];
 171                 r_nn[2][i] = r2;
 172                 nn[2][i]   = j;
 173             }
 174             else if (r2 < r_nn[3][i])
 175             {
 176                 r_nn[3][i] = r2;
 177                 nn[3][i]   = j;
 178             }
 179         }
 180
 181
 182         /* calculate mean distance between nearest neighbours */
 183         rmean = 0;
 184         for (j = 0; (j < 4); j++)
 185         {
 186             r_nn[j][i] = std::sqrt(r_nn[j][i]);
 187             rmean += r_nn[j][i];
 188         }
 189         rmean /= 4;
 190
 191         n        = 0;
 192         sgmol[i] = 0.0;
 193         skmol[i] = 0.0;
 194
 195         /* Chau1998a eqn 3 */
 196         /* angular part tetrahedrality order parameter per atom */
 197         for (j = 0; (j < 3); j++)
 198         {
 199             for (k = j + 1; (k < 4); k++)
 200             {
 201                 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
 202                 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
 203
 204                 unitv(rk, urk);
 205                 unitv(rj, urj);
 206
 207                 cost  = iprod(urk, urj) + onethird;
 208                 cost2 = cost * cost;
 209
 210                 /* sgmol[i] += 3*cost2/32;  */
 211                 sgmol[i] += cost2;
 212
 213                 /* determine distribution */
 214                 ibin = static_cast<int>(static_cast<real>(nsgbin) * cost2);
 215                 if (ibin < nsgbin)
 216                 {
 217                     sgbin[ibin]++;
 218                 }
 219                 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
 220                 l++;
 221                 n++;
 222             }
 223         }
 224
 225         /* normalize sgmol between 0.0 and 1.0 */
 226         sgmol[i] = 3 * sgmol[i] / 32;
 227         *sgmean += sgmol[i];
 228
 229         /* distance part tetrahedrality order parameter per atom */
 230         rmean2 = 4 * 3 * rmean * rmean;
 231         for (j = 0; (j < 4); j++)
 232         {
 233             skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
 234             /*      printf("%d %f (%f %f %f %f) \n",
 235                 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
 236              */
 237         }
 238
 239         *skmean += skmol[i];
 240
 241         /* Compute sliced stuff */
 242         sl_index = static_cast<int>(std::round((1 + x[i][slice_dim] / box[slice_dim][slice_dim])
 243                                                * static_cast<real>(nslice)))
 244                    % nslice;
 245         sgslice[sl_index] += sgmol[i];
 246         skslice[sl_index] += skmol[i];
 247         sl_count[sl_index]++;
 248     } /* loop over entries in index file */
 249
 250     *sgmean /= static_cast<real>(maxidx);
 251     *skmean /= static_cast<real>(maxidx);
 252
 253     for (i = 0; (i < nslice); i++)
 254     {
 255         if (sl_count[i] > 0)
 256         {
 257             sgslice[i] /= sl_count[i];
 258             skslice[i] /= sl_count[i];
 259         }
 260     }
 261     sfree(sl_count);
 262     sfree(sgbin);
 263     sfree(sgmol);
 264     sfree(skmol);
 265     for (i = 0; (i < 4); i++)
 266     {
 267         sfree(r_nn[i]);
 268         sfree(nn[i]);
 269     }
 270 }
 271
 272
 273 static void calc_tetra_order_parm(const char*             fnNDX,
 274                                   const char*             fnTPS,
 275                                   const char*             fnTRX,
 276                                   const char*             sgfn,
 277                                   const char*             skfn,
 278                                   int                     nslice,
 279                                   int                     slice_dim,
 280                                   const char*             sgslfn,
 281                                   const char*             skslfn,
 282                                   const gmx_output_env_t* oenv)
 283 {
 284     FILE *       fpsg = nullptr, *fpsk = nullptr;
 285     t_topology   top;
 286     PbcType      pbcType;
 287     t_trxstatus* status;
 288     int          natoms;
 289     real         t;
 290     rvec *       xtop, *x;
 291     matrix       box;
 292     real         sg, sk;
 293     int**        index;
 294     char**       grpname;
 295     int          i, *isize, ng, nframes;
 296     real *       sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
 297     gmx_rmpbc_t  gpbc = nullptr;
 298
 299
 300     read_tps_conf(fnTPS, &top, &pbcType, &xtop, nullptr, box, FALSE);
 301
 302     snew(sg_slice, nslice);
 303     snew(sk_slice, nslice);
 304     snew(sg_slice_tot, nslice);
 305     snew(sk_slice_tot, nslice);
 306     ng = 1;
 307     /* get index groups */
 308     printf("Select the group that contains the atoms you want to use for the tetrahedrality order "
 309            "parameter calculation:\n");
 310     snew(grpname, ng);
 311     snew(index, ng);
 312     snew(isize, ng);
 313     get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
 314
 315     /* Analyze trajectory */
 316     natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
 317     if (natoms > top.atoms.nr)
 318     {
 319         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)", top.atoms.nr, natoms);
 320     }
 321     check_index(nullptr, ng, index[0], nullptr, natoms);
 322
 323     fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N", oenv);
 324     fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N", oenv);
 325
 326     /* loop over frames */
 327     gpbc    = gmx_rmpbc_init(&top.idef, pbcType, natoms);
 328     nframes = 0;
 329     do
 330     {
 331         find_nearest_neighbours(
 332                 pbcType, natoms, box, x, isize[0], index[0], &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
 333         for (i = 0; (i < nslice); i++)
 334         {
 335             sg_slice_tot[i] += sg_slice[i];
 336             sk_slice_tot[i] += sk_slice[i];
 337         }
 338         fprintf(fpsg, "%f %f\n", t, sg);
 339         fprintf(fpsk, "%f %f\n", t, sk);
 340         nframes++;
 341     } while (read_next_x(oenv, status, &t, x, box));
 342     close_trx(status);
 343     gmx_rmpbc_done(gpbc);
 344
 345     sfree(grpname);
 346     sfree(index);
 347     sfree(isize);
 348
 349     xvgrclose(fpsg);
 350     xvgrclose(fpsk);
 351
 352     fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N", oenv);
 353     fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N", oenv);
 354     for (i = 0; (i < nslice); i++)
 355     {
 356         fprintf(fpsg,
 357                 "%10g  %10g\n",
 358                 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
 359                 sg_slice_tot[i] / static_cast<real>(nframes));
 360         fprintf(fpsk,
 361                 "%10g  %10g\n",
 362                 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
 363                 sk_slice_tot[i] / static_cast<real>(nframes));
 364     }
 365     xvgrclose(fpsg);
 366     xvgrclose(fpsk);
 367 }
 368
 369
 370 /* Print name of first atom in all groups in index file */
 371 static void print_types(const int index[], int a[], int ngrps, char* groups[], const t_topology* top)
 372 {
 373     int i;
 374
 375     fprintf(stderr, "Using following groups: \n");
 376     for (i = 0; i < ngrps; i++)
 377     {
 378         fprintf(stderr,
 379                 "Groupname: %s First atomname: %s First atomnr %d\n",
 380                 groups[i],
 381                 *(top->atoms.atomname[a[index[i]]]),
 382                 a[index[i]]);
 383     }
 384     fprintf(stderr, "\n");
 385 }
 386
 387 static void check_length(real length, int a, int b)
 388 {
 389     if (length > 0.3)
 390     {
 391         fprintf(stderr,
 392                 "WARNING: distance between atoms %d and "
 393                 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
 394                 a,
 395                 b,
 396                 length);
 397     }
 398 }
 399
 400 static void calc_order(const char*             fn,
 401                        const int*              index,
 402                        int*                    a,
 403                        rvec**                  order,
 404                        real***                 slOrder,
 405                        real*                   slWidth,
 406                        int                     nslices,
 407                        gmx_bool                bSliced,
 408                        gmx_bool                bUnsat,
 409                        const t_topology*       top,
 410                        PbcType                 pbcType,
 411                        int                     ngrps,
 412                        int                     axis,
 413                        gmx_bool                permolecule,
 414                        gmx_bool                radial,
 415                        gmx_bool                distcalc,
 416                        const char*             radfn,
 417                        real***                 distvals,
 418                        const gmx_output_env_t* oenv)
 419 {
 420     /* if permolecule = TRUE, order parameters will be calculed per molecule
 421      * and stored in slOrder with #slices = # molecules */
 422     rvec *x0,         /* coordinates with pbc                           */
 423             *x1;      /* coordinates without pbc                        */
 424     matrix       box; /* box (3x3)                                      */
 425     t_trxstatus* status;
 426     rvec         cossum,                   /* sum of vector angles for three axes            */
 427             Sx, Sy, Sz,                    /* the three molecular axes                       */
 428             tmp1, tmp2,                    /* temp. rvecs for calculating dot products       */
 429             frameorder;                    /* order parameters for one frame                 */
 430     real* slFrameorder;                    /* order parameter for one frame, per slice      */
 431     real  length,                          /* total distance between two atoms               */
 432             t,                             /* time from trajectory                           */
 433             z_ave, z1, z2;                 /* average z, used to det. which slice atom is in */
 434     int natoms,                            /* nr. atoms in trj                               */
 435             nr_tails,                      /* nr tails, to check if index file is correct    */
 436             size = 0,                      /* nr. of atoms in group. same as nr_tails        */
 437             i, j, m, k, teller = 0, slice; /* current slice number                           */
 438     real        nr_frames = 0;
 439     int*        slCount;      /* nr. of atoms in one slice                      */
 440     real        sdbangle = 0; /* sum of these angles                            */
 441     gmx_bool    use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
 442     rvec        direction, com;
 443     int         comsize, distsize;
 444     int *       comidx = nullptr, *distidx = nullptr;
 445     char*       grpname = nullptr;
 446     t_pbc       pbc;
 447     real        arcdist, tmpdist;
 448     gmx_rmpbc_t gpbc = nullptr;
 449
 450     /* PBC added for center-of-mass vector*/
 451     /* Initiate the pbc structure */
 452     std::memset(&pbc, 0, sizeof(pbc));
 453
 454     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 455     {
 456         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 457     }
 458
 459     nr_tails = index[1] - index[0];
 460     fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
 461     /* take first group as standard. Not rocksolid, but might catch error in index*/
 462
 463     if (permolecule)
 464     {
 465         nslices = nr_tails;
 466         bSliced = FALSE; /*force slices off */
 467         fprintf(stderr, "Calculating order parameters for each of %d molecules\n", nslices);
 468     }
 469
 470     if (radial)
 471     {
 472         use_unitvector = TRUE;
 473         fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
 474         get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
 475     }
 476     if (distcalc)
 477     {
 478         if (grpname != nullptr)
 479         {
 480             sfree(grpname);
 481         }
 482         fprintf(stderr, "Select an index group to use as distance reference\n");
 483         get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
 484         bSliced = FALSE; /*force slices off*/
 485     }
 486
 487     if (use_unitvector && bSliced)
 488     {
 489         fprintf(stderr,
 490                 "Warning:  slicing and specified unit vectors are not currently compatible\n");
 491     }
 492
 493     snew(slCount, nslices);
 494     snew(*slOrder, nslices);
 495     for (i = 0; i < nslices; i++)
 496     {
 497         snew((*slOrder)[i], ngrps);
 498     }
 499     if (distcalc)
 500     {
 501         snew(*distvals, nslices);
 502         for (i = 0; i < nslices; i++)
 503         {
 504             snew((*distvals)[i], ngrps);
 505         }
 506     }
 507     snew(*order, ngrps);
 508     snew(slFrameorder, nslices);
 509     snew(x1, natoms);
 510
 511     if (bSliced)
 512     {
 513         *slWidth = box[axis][axis] / static_cast<real>(nslices);
 514         fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", nslices, *slWidth);
 515     }
 516
 517
 518 #if 0
 519     nr_tails = index[1] - index[0];
 520     fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
 521     /* take first group as standard. Not rocksolid, but might catch error
 522        in index*/
 523 #endif
 524
 525     teller = 0;
 526
 527     gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
 528     /*********** Start processing trajectory ***********/
 529     do
 530     {
 531         if (bSliced)
 532         {
 533             *slWidth = box[axis][axis] / static_cast<real>(nslices);
 534         }
 535         teller++;
 536
 537         set_pbc(&pbc, pbcType, box);
 538         gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
 539
 540         /* Now loop over all groups. There are ngrps groups, the order parameter can
 541            be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
 542            atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
 543            course, in this case index[i+1] -index[i] has to be the same for all
 544            groups, namely the number of tails. i just runs over all atoms in a tail,
 545            so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
 546          */
 547
 548
 549         if (radial)
 550         {
 551             /*center-of-mass determination*/
 552             com[XX] = 0.0;
 553             com[YY] = 0.0;
 554             com[ZZ] = 0.0;
 555             for (j = 0; j < comsize; j++)
 556             {
 557                 rvec_inc(com, x1[comidx[j]]);
 558             }
 559             svmul(1.0 / comsize, com, com);
 560         }
 561         rvec displacementFromReference;
 562         if (distcalc)
 563         {
 564             rvec dref = { 0.0, 0.0, 0.0 };
 565             for (j = 0; j < distsize; j++)
 566             {
 567                 rvec_inc(dref, x1[distidx[j]]);
 568             }
 569             svmul(1.0 / distsize, dref, dref);
 570             if (radial)
 571             {
 572                 pbc_dx(&pbc, dref, com, displacementFromReference);
 573                 unitv(displacementFromReference, displacementFromReference);
 574             }
 575         }
 576
 577         for (i = 1; i < ngrps - 1; i++)
 578         {
 579             clear_rvec(frameorder);
 580
 581             size = index[i + 1] - index[i];
 582             if (size != nr_tails)
 583             {
 584                 gmx_fatal(FARGS,
 585                           "grp %d does not have same number of"
 586                           " elements as grp 1\n",
 587                           i);
 588             }
 589
 590             for (j = 0; j < size; j++)
 591             {
 592                 if (radial)
 593                 /*create unit vector*/
 594                 {
 595                     pbc_dx(&pbc, x1[a[index[i] + j]], com, direction);
 596                     unitv(direction, direction);
 597                     /*DEBUG*/
 598                     /*if (j==0)
 599                         fprintf(stderr,"X %f %f %f\tcom %f %f %f\tdirection %f %f %f\n",x1[a[index[i]+j]][0],x1[a[index[i]+j]][1],x1[a[index[i]+j]][2],com[0],com[1],com[2],
 600                             direction[0],direction[1],direction[2]);*/
 601                 }
 602
 603                 if (bUnsat)
 604                 {
 605                     rvec dist;
 606                     /* Using convention for unsaturated carbons */
 607                     /* first get Sz, the vector from Cn to Cn+1 */
 608                     rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], dist);
 609                     length = norm(dist);
 610                     check_length(length, a[index[i] + j], a[index[i + 1] + j]);
 611                     svmul(1.0 / length, dist, Sz);
 612
 613                     /* this is actually the cosine of the angle between the double bond
 614                        and axis, because Sz is normalized and the two other components of
 615                        the axis on the bilayer are zero */
 616                     if (use_unitvector)
 617                     {
 618                         sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
 619                     }
 620                     else
 621                     {
 622                         sdbangle += std::acos(Sz[axis]);
 623                     }
 624                 }
 625                 else
 626                 {
 627                     rvec dist;
 628                     /* get vector dist(Cn-1,Cn+1) for tail atoms */
 629                     rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i - 1] + j]], dist);
 630                     length = norm(dist); /* determine distance between two atoms */
 631                     check_length(length, a[index[i - 1] + j], a[index[i + 1] + j]);
 632
 633                     svmul(1.0 / length, dist, Sz);
 634                     /* Sz is now the molecular axis Sz, normalized and all that */
 635                 }
 636
 637                 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
 638                    we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
 639                 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], tmp1);
 640                 rvec_sub(x1[a[index[i - 1] + j]], x1[a[index[i] + j]], tmp2);
 641                 cprod(tmp1, tmp2, Sx);
 642                 svmul(1.0 / norm(Sx), Sx, Sx);
 643
 644                 /* now we can get Sy from the outer product of Sx and Sz   */
 645                 cprod(Sz, Sx, Sy);
 646                 svmul(1.0 / norm(Sy), Sy, Sy);
 647
 648                 /* the square of cosine of the angle between dist and the axis.
 649                    Using the innerproduct, but two of the three elements are zero
 650                    Determine the sum of the orderparameter of all atoms in group
 651                  */
 652                 if (use_unitvector)
 653                 {
 654                     cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
 655                     cossum[YY] = gmx::square(iprod(Sy, direction));
 656                     cossum[ZZ] = gmx::square(iprod(Sz, direction));
 657                 }
 658                 else
 659                 {
 660                     cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
 661                     cossum[YY] = gmx::square(Sy[axis]);
 662                     cossum[ZZ] = gmx::square(Sz[axis]);
 663                 }
 664
 665                 for (m = 0; m < DIM; m++)
 666                 {
 667                     frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
 668                 }
 669
 670                 if (bSliced)
 671                 {
 672                     /* get average coordinate in box length for slicing,
 673                        determine which slice atom is in, increase count for that
 674                        slice. slFrameorder and slOrder are reals, not
 675                        rvecs. Only the component [axis] of the order tensor is
 676                        kept, until I find it necessary to know the others too
 677                      */
 678
 679                     z1    = x1[a[index[i - 1] + j]][axis];
 680                     z2    = x1[a[index[i + 1] + j]][axis];
 681                     z_ave = 0.5 * (z1 + z2);
 682                     slice = static_cast<int>((static_cast<real>(nslices) * z_ave) / box[axis][axis]);
 683                     while (slice < 0)
 684                     {
 685                         slice += static_cast<real>(nslices);
 686                     }
 687                     slice = slice % nslices;
 688
 689                     slCount[slice]++; /* determine slice, increase count */
 690
 691                     slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
 692                 }
 693                 else if (permolecule)
 694                 {
 695                     /*  store per-molecule order parameter
 696                      *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 *
 697                      * iprod(cossum,direction) - 1); following is for Scd order: */
 698                     (*slOrder)[j][i] += -1
 699                                         * (1.0 / 3.0 * (3 * cossum[XX] - 1)
 700                                            + 1.0 / 3.0 * 0.5 * (3.0 * cossum[YY] - 1));
 701                 }
 702                 if (distcalc)
 703                 {
 704                     if (radial)
 705                     {
 706                         /* bin order parameter by arc distance from reference group*/
 707                         arcdist = gmx_angle(displacementFromReference, direction);
 708                         (*distvals)[j][i] += arcdist;
 709                     }
 710                     else if (i == 1)
 711                     {
 712                         /* Want minimum lateral distance to first group calculated */
 713                         tmpdist = trace(box); /* should be max value */
 714                         for (k = 0; k < distsize; k++)
 715                         {
 716                             rvec displacement;
 717                             pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i] + j]], displacement);
 718                             /* at the moment, just remove displacement[axis] */
 719                             displacement[axis] = 0;
 720                             tmpdist            = std::min(tmpdist, norm2(displacement));
 721                         }
 722                         // fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
 723                         (*distvals)[j][i] += std::sqrt(tmpdist);
 724                     }
 725                 }
 726             } /* end loop j, over all atoms in group */
 727
 728             for (m = 0; m < DIM; m++)
 729             {
 730                 (*order)[i][m] += (frameorder[m] / static_cast<real>(size));
 731             }
 732
 733             if (!permolecule)
 734             { /*Skip following if doing per-molecule*/
 735                 for (k = 0; k < nslices; k++)
 736                 {
 737                     if (slCount[k]) /* if no elements, nothing has to be added */
 738                     {
 739                         (*slOrder)[k][i] += slFrameorder[k] / static_cast<real>(slCount[k]);
 740                         slFrameorder[k] = 0;
 741                         slCount[k]      = 0;
 742                     }
 743                 }
 744             } /* end loop i, over all groups in indexfile */
 745         }
 746         nr_frames++;
 747
 748     } while (read_next_x(oenv, status, &t, x0, box));
 749     /*********** done with status file **********/
 750
 751     fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
 752     gmx_rmpbc_done(gpbc);
 753
 754     /* average over frames */
 755     for (i = 1; i < ngrps - 1; i++)
 756     {
 757         svmul(1.0 / nr_frames, (*order)[i], (*order)[i]);
 758         fprintf(stderr,
 759                 "Atom %d Tensor: x=%g , y=%g, z=%g\n",
 760                 i,
 761                 (*order)[i][XX],
 762                 (*order)[i][YY],
 763                 (*order)[i][ZZ]);
 764         if (bSliced || permolecule)
 765         {
 766             for (k = 0; k < nslices; k++)
 767             {
 768                 (*slOrder)[k][i] /= nr_frames;
 769             }
 770         }
 771         if (distcalc)
 772         {
 773             for (k = 0; k < nslices; k++)
 774             {
 775                 (*distvals)[k][i] /= nr_frames;
 776             }
 777         }
 778     }
 779
 780     if (bUnsat)
 781     {
 782         fprintf(stderr,
 783                 "Average angle between double bond and normal: %f\n",
 784                 180 * sdbangle / (nr_frames * static_cast<real>(size) * M_PI));
 785     }
 786
 787     sfree(x0); /* free memory used by coordinate arrays */
 788     sfree(x1);
 789     if (comidx != nullptr)
 790     {
 791         sfree(comidx);
 792     }
 793     if (distidx != nullptr)
 794     {
 795         sfree(distidx);
 796     }
 797     if (grpname != nullptr)
 798     {
 799         sfree(grpname);
 800     }
 801 }
 802
 803
 804 static void order_plot(rvec                    order[],
 805                        real*                   slOrder[],
 806                        const char*             afile,
 807                        const char*             bfile,
 808                        const char*             cfile,
 809                        int                     ngrps,
 810                        int                     nslices,
 811                        real                    slWidth,
 812                        gmx_bool                bSzonly,
 813                        gmx_bool                permolecule,
 814                        real**                  distvals,
 815                        const gmx_output_env_t* oenv)
 816 {
 817     FILE *ord, *slOrd; /* xvgr files with order parameters  */
 818     int   atom, slice; /* atom corresponding to order para.*/
 819     char  buf[256];    /* for xvgr title */
 820     real  S;           /* order parameter averaged over all atoms */
 821
 822     if (permolecule)
 823     {
 824         sprintf(buf, "Scd order parameters");
 825         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 826         sprintf(buf, "Orderparameters per atom per slice");
 827         slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
 828         for (atom = 1; atom < ngrps - 1; atom++)
 829         {
 830             fprintf(ord,
 831                     "%12d   %12g\n",
 832                     atom,
 833                     -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
 834         }
 835
 836         for (slice = 0; slice < nslices; slice++)
 837         {
 838             fprintf(slOrd, "%12d\t", slice);
 839             if (distvals)
 840             {
 841                 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
 842             }
 843             for (atom = 1; atom < ngrps - 1; atom++)
 844             {
 845                 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
 846             }
 847             fprintf(slOrd, "\n");
 848         }
 849     }
 850     else if (bSzonly)
 851     {
 852         sprintf(buf, "Orderparameters Sz per atom");
 853         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 854         fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
 855
 856         sprintf(buf, "Orderparameters per atom per slice");
 857         slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
 858
 859         for (atom = 1; atom < ngrps - 1; atom++)
 860         {
 861             fprintf(ord, "%12d       %12g\n", atom, order[atom][ZZ]);
 862         }
 863
 864         for (slice = 0; slice < nslices; slice++)
 865         {
 866             S = 0;
 867             for (atom = 1; atom < ngrps - 1; atom++)
 868             {
 869                 S += slOrder[slice][atom];
 870             }
 871             fprintf(slOrd, "%12g     %12g\n", static_cast<real>(slice) * slWidth, S / static_cast<real>(atom));
 872         }
 873     }
 874     else
 875     {
 876         sprintf(buf, "Order tensor diagonal elements");
 877         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 878         sprintf(buf, "Deuterium order parameters");
 879         slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
 880
 881         for (atom = 1; atom < ngrps - 1; atom++)
 882         {
 883             fprintf(ord, "%12d   %12g   %12g   %12g\n", atom, order[atom][XX], order[atom][YY], order[atom][ZZ]);
 884             fprintf(slOrd,
 885                     "%12d   %12g\n",
 886                     atom,
 887                     -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
 888         }
 889     }
 890     xvgrclose(ord);
 891     xvgrclose(slOrd);
 892 }
 893
 894 static void write_bfactors(t_filenm*         fnm,
 895                            int               nfile,
 896                            const int*        index,
 897                            const int*        a,
 898                            int               nslices,
 899                            int               ngrps,
 900                            real**            order,
 901                            const t_topology* top,
 902                            real**            distvals,
 903                            gmx_output_env_t* oenv)
 904 {
 905     /*function to write order parameters as B factors in PDB file using
 906           first frame of trajectory*/
 907     t_trxstatus* status;
 908     t_trxframe   fr, frout;
 909     t_atoms      useatoms;
 910     int          i, j, ctr, nout;
 911
 912     ngrps -= 2; /*we don't have an order parameter for the first or
 913                       last atom in each chain*/
 914     nout = nslices * ngrps;
 915     read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
 916
 917     close_trx(status);
 918     frout        = fr;
 919     frout.natoms = nout;
 920     frout.bF     = FALSE;
 921     frout.bV     = FALSE;
 922     frout.x      = nullptr;
 923     snew(frout.x, nout);
 924
 925     init_t_atoms(&useatoms, nout, TRUE);
 926     useatoms.nr = nout;
 927
 928     /*initialize PDBinfo*/
 929     for (i = 0; i < useatoms.nr; ++i)
 930     {
 931         useatoms.pdbinfo[i].type         = 0;
 932         useatoms.pdbinfo[i].occup        = 0.0;
 933         useatoms.pdbinfo[i].bfac         = 0.0;
 934         useatoms.pdbinfo[i].bAnisotropic = FALSE;
 935     }
 936
 937     for (j = 0, ctr = 0; j < nslices; j++)
 938     {
 939         for (i = 0; i < ngrps; i++, ctr++)
 940         {
 941             /*iterate along each chain*/
 942             useatoms.pdbinfo[ctr].bfac = order[j][i + 1];
 943             if (distvals)
 944             {
 945                 useatoms.pdbinfo[ctr].occup = distvals[j][i + 1];
 946             }
 947             copy_rvec(fr.x[a[index[i + 1] + j]], frout.x[ctr]);
 948             useatoms.atomname[ctr] = top->atoms.atomname[a[index[i + 1] + j]];
 949             useatoms.atom[ctr]     = top->atoms.atom[a[index[i + 1] + j]];
 950             useatoms.nres          = std::max(useatoms.nres, useatoms.atom[ctr].resind + 1);
 951             useatoms.resinfo[useatoms.atom[ctr].resind] =
 952                     top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
 953         }
 954     }
 955
 956     write_sto_conf(
 957             opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, nullptr, frout.pbcType, frout.box);
 958
 959     sfree(frout.x);
 960     done_atom(&useatoms);
 961 }
 962
 963 int gmx_order(int argc, char* argv[])
 964 {
 965     const char* desc[] = {
 966         "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
 967         "vector i-1, i+1 is used together with an axis. ",
 968         "The index file should contain only the groups to be used for calculations,",
 969         "with each group of equivalent carbons along the relevant acyl chain in its own",
 970         "group. There should not be any generic groups (like System, Protein) in the index",
 971         "file to avoid confusing the program (this is not relevant to tetrahedral order",
 972         "parameters however, which only work for water anyway).[PAR]",
 973         "[THISMODULE] can also give all",
 974         "diagonal elements of the order tensor and even calculate the deuterium",
 975         "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
 976         "order tensor component (specified by the [TT]-d[tt] option) is given and the",
 977         "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
 978         "selected, all diagonal elements and the deuterium order parameter is",
 979         "given.[PAR]",
 980         "The tetrahedrality order parameters can be determined",
 981         "around an atom. Both angle an distance order parameters are calculated. See",
 982         "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
 983         "for more details."
 984     };
 985
 986     static int         nslices       = 1;     /* nr of slices defined       */
 987     static gmx_bool    bSzonly       = FALSE; /* True if only Sz is wanted  */
 988     static gmx_bool    bUnsat        = FALSE; /* True if carbons are unsat. */
 989     static const char* normal_axis[] = { nullptr, "z", "x", "y", nullptr };
 990     static gmx_bool    permolecule   = FALSE; /*compute on a per-molecule basis */
 991     static gmx_bool    radial        = FALSE; /*compute a radial membrane normal */
 992     static gmx_bool    distcalc      = FALSE; /*calculate distance from a reference group */
 993     t_pargs            pa[]          = {
 994         { "-d", FALSE, etENUM, { normal_axis }, "Direction of the normal on the membrane" },
 995         { "-sl",
 996           FALSE,
 997           etINT,
 998           { &nslices },
 999           "Calculate order parameter as function of box length, dividing the box"
1000           " into this number of slices." },
1001         { "-szonly",
1002           FALSE,
1003           etBOOL,
1004           { &bSzonly },
1005           "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
1006         { "-unsat",
1007           FALSE,
1008           etBOOL,
1009           { &bUnsat },
1010           "Calculate order parameters for unsaturated carbons. Note that this can"
1011           "not be mixed with normal order parameters." },
1012         { "-permolecule",
1013           FALSE,
1014           etBOOL,
1015           { &permolecule },
1016           "Compute per-molecule Scd order parameters" },
1017         { "-radial", FALSE, etBOOL, { &radial }, "Compute a radial membrane normal" },
1018         { "-calcdist", FALSE, etBOOL, { &distcalc }, "Compute distance from a reference" },
1019     };
1020
1021     rvec*  order;         /* order par. for each atom   */
1022     real** slOrder;       /* same, per slice            */
1023     real   slWidth = 0.0; /* width of a slice           */
1024     char** grpname;       /* groupnames                 */
1025     int    ngrps,         /* nr. of groups              */
1026             i, axis = 0;  /* normal axis                */
1027     t_topology* top;      /* topology                   */
1028     PbcType     pbcType;  /* type of periodic boundary conditions */
1029     int *       index,    /* indices for a              */
1030             *a;           /* atom numbers in each group */
1031     t_blocka* block;      /* data from index file       */
1032     t_filenm  fnm[] = {
1033         /* files for g_order    */
1034         { efTRX, "-f", nullptr, ffREAD },     /* trajectory file              */
1035         { efNDX, "-n", nullptr, ffREAD },     /* index file           */
1036         { efNDX, "-nr", nullptr, ffOPTRD },   /* index for radial axis calculation */
1037         { efTPR, nullptr, nullptr, ffREAD },  /* topology file                */
1038         { efXVG, "-o", "order", ffWRITE },    /* xvgr output file     */
1039         { efXVG, "-od", "deuter", ffWRITE },  /* xvgr output file           */
1040         { efPDB, "-ob", nullptr, ffOPTWR },   /* write Scd as B factors to PDB if permolecule   */
1041         { efXVG, "-os", "sliced", ffWRITE },  /* xvgr output file           */
1042         { efXVG, "-Sg", "sg-ang", ffOPTWR },  /* xvgr output file           */
1043         { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file           */
1044         { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR },  /* xvgr output file           */
1045         { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file           */
1046     };
1047     gmx_bool bSliced = FALSE; /* True if box is sliced      */
1048 #define NFILE asize(fnm)
1049     real**            distvals = nullptr;
1050     const char *      sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
1051     gmx_output_env_t* oenv;
1052
1053     if (!parse_common_args(
1054                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1055     {
1056         return 0;
1057     }
1058     if (nslices < 1)
1059     {
1060         gmx_fatal(FARGS, "Can not have nslices < 1");
1061     }
1062     sgfnm  = opt2fn_null("-Sg", NFILE, fnm);
1063     skfnm  = opt2fn_null("-Sk", NFILE, fnm);
1064     ndxfnm = opt2fn_null("-n", NFILE, fnm);
1065     tpsfnm = ftp2fn(efTPR, NFILE, fnm);
1066     trxfnm = ftp2fn(efTRX, NFILE, fnm);
1067
1068     /* Calculate axis */
1069     GMX_RELEASE_ASSERT(normal_axis[0] != nullptr, "Options inconsistency; normal_axis[0] is NULL");
1070     if (std::strcmp(normal_axis[0], "x") == 0)
1071     {
1072         axis = XX;
1073     }
1074     else if (std::strcmp(normal_axis[0], "y") == 0)
1075     {
1076         axis = YY;
1077     }
1078     else if (std::strcmp(normal_axis[0], "z") == 0)
1079     {
1080         axis = ZZ;
1081     }
1082     else
1083     {
1084         gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1085     }
1086
1087     switch (axis)
1088     {
1089         case 0: fprintf(stderr, "Taking x axis as normal to the membrane\n"); break;
1090         case 1: fprintf(stderr, "Taking y axis as normal to the membrane\n"); break;
1091         case 2: fprintf(stderr, "Taking z axis as normal to the membrane\n"); break;
1092     }
1093
1094     /* tetraheder order parameter */
1095     if (skfnm || sgfnm)
1096     {
1097         /* If either of theoptions is set we compute both */
1098         sgfnm = opt2fn("-Sg", NFILE, fnm);
1099         skfnm = opt2fn("-Sk", NFILE, fnm);
1100         calc_tetra_order_parm(ndxfnm,
1101                               tpsfnm,
1102                               trxfnm,
1103                               sgfnm,
1104                               skfnm,
1105                               nslices,
1106                               axis,
1107                               opt2fn("-Sgsl", NFILE, fnm),
1108                               opt2fn("-Sksl", NFILE, fnm),
1109                               oenv);
1110         /* view xvgr files */
1111         do_view(oenv, opt2fn("-Sg", NFILE, fnm), nullptr);
1112         do_view(oenv, opt2fn("-Sk", NFILE, fnm), nullptr);
1113         if (nslices > 1)
1114         {
1115             do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), nullptr);
1116             do_view(oenv, opt2fn("-Sksl", NFILE, fnm), nullptr);
1117         }
1118     }
1119     else
1120     {
1121         /* tail order parameter */
1122
1123         if (nslices > 1)
1124         {
1125             bSliced = TRUE;
1126             fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1127         }
1128
1129         if (bSzonly)
1130         {
1131             fprintf(stderr, "Only calculating Sz\n");
1132         }
1133         if (bUnsat)
1134         {
1135             fprintf(stderr, "Taking carbons as unsaturated!\n");
1136         }
1137
1138         top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
1139
1140         block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1141         index = block->index; /* get indices from t_block block */
1142         a     = block->a;     /* see block.h                    */
1143         ngrps = block->nr;
1144
1145         if (permolecule)
1146         {
1147             nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1148             fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1149         }
1150
1151         if (radial)
1152         {
1153             fprintf(stderr, "Calculating radial distances\n");
1154             if (!permolecule)
1155             {
1156                 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1157             }
1158         }
1159
1160         /* show atomtypes, to check if index file is correct */
1161         print_types(index, a, ngrps, grpname, top);
1162
1163         calc_order(ftp2fn(efTRX, NFILE, fnm),
1164                    index,
1165                    a,
1166                    &order,
1167                    &slOrder,
1168                    &slWidth,
1169                    nslices,
1170                    bSliced,
1171                    bUnsat,
1172                    top,
1173                    pbcType,
1174                    ngrps,
1175                    axis,
1176                    permolecule,
1177                    radial,
1178                    distcalc,
1179                    opt2fn_null("-nr", NFILE, fnm),
1180                    &distvals,
1181                    oenv);
1182
1183         if (radial)
1184         {
1185             ngrps--; /*don't print the last group--was used for
1186                                center-of-mass determination*/
1187         }
1188         order_plot(order,
1189                    slOrder,
1190                    opt2fn("-o", NFILE, fnm),
1191                    opt2fn("-os", NFILE, fnm),
1192                    opt2fn("-od", NFILE, fnm),
1193                    ngrps,
1194                    nslices,
1195                    slWidth,
1196                    bSzonly,
1197                    permolecule,
1198                    distvals,
1199                    oenv);
1200
1201         if (opt2bSet("-ob", NFILE, fnm))
1202         {
1203             if (!permolecule)
1204             {
1205                 fprintf(stderr,
1206                         "Won't write B-factors with averaged order parameters; use -permolecule\n");
1207             }
1208             else
1209             {
1210                 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1211             }
1212         }
1213
1214
1215         do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);  /* view xvgr file */
1216         do_view(oenv, opt2fn("-os", NFILE, fnm), nullptr); /* view xvgr file */
1217         do_view(oenv, opt2fn("-od", NFILE, fnm), nullptr); /* view xvgr file */
1218     }
1219
1220     if (distvals != nullptr)
1221     {
1222         for (i = 0; i < nslices; ++i)
1223         {
1224             sfree(distvals[i]);
1225         }
1226         sfree(distvals);
1227     }
1228
1229     return 0;
1230 }