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