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(int         ePBC,
  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, ePBC, 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     int          ePBC;
 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, &ePBC, &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, ePBC, natoms);
 328     nframes = 0;
 329     do
 330     {
 331         find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0], &sg, &sk, nslice,
 332                                 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, "%10g  %10g\n", (i + 0.5) * box[slice_dim][slice_dim] / nslice,
 357                 sg_slice_tot[i] / static_cast<real>(nframes));
 358         fprintf(fpsk, "%10g  %10g\n", (i + 0.5) * box[slice_dim][slice_dim] / nslice,
 359                 sk_slice_tot[i] / static_cast<real>(nframes));
 360     }
 361     xvgrclose(fpsg);
 362     xvgrclose(fpsk);
 363 }
 364
 365
 366 /* Print name of first atom in all groups in index file */
 367 static void print_types(const int index[], int a[], int ngrps, char* groups[], const t_topology* top)
 368 {
 369     int i;
 370
 371     fprintf(stderr, "Using following groups: \n");
 372     for (i = 0; i < ngrps; i++)
 373     {
 374         fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %d\n", groups[i],
 375                 *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
 376     }
 377     fprintf(stderr, "\n");
 378 }
 379
 380 static void check_length(real length, int a, int b)
 381 {
 382     if (length > 0.3)
 383     {
 384         fprintf(stderr,
 385                 "WARNING: distance between atoms %d and "
 386                 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
 387                 a, b, length);
 388     }
 389 }
 390
 391 static void calc_order(const char*             fn,
 392                        const int*              index,
 393                        int*                    a,
 394                        rvec**                  order,
 395                        real***                 slOrder,
 396                        real*                   slWidth,
 397                        int                     nslices,
 398                        gmx_bool                bSliced,
 399                        gmx_bool                bUnsat,
 400                        const t_topology*       top,
 401                        int                     ePBC,
 402                        int                     ngrps,
 403                        int                     axis,
 404                        gmx_bool                permolecule,
 405                        gmx_bool                radial,
 406                        gmx_bool                distcalc,
 407                        const char*             radfn,
 408                        real***                 distvals,
 409                        const gmx_output_env_t* oenv)
 410 {
 411     /* if permolecule = TRUE, order parameters will be calculed per molecule
 412      * and stored in slOrder with #slices = # molecules */
 413     rvec *x0,         /* coordinates with pbc                           */
 414             *x1;      /* coordinates without pbc                        */
 415     matrix       box; /* box (3x3)                                      */
 416     t_trxstatus* status;
 417     rvec         cossum,                   /* sum of vector angles for three axes            */
 418             Sx, Sy, Sz,                    /* the three molecular axes                       */
 419             tmp1, tmp2,                    /* temp. rvecs for calculating dot products       */
 420             frameorder;                    /* order parameters for one frame                 */
 421     real* slFrameorder;                    /* order parameter for one frame, per slice      */
 422     real  length,                          /* total distance between two atoms               */
 423             t,                             /* time from trajectory                           */
 424             z_ave, z1, z2;                 /* average z, used to det. which slice atom is in */
 425     int natoms,                            /* nr. atoms in trj                               */
 426             nr_tails,                      /* nr tails, to check if index file is correct    */
 427             size = 0,                      /* nr. of atoms in group. same as nr_tails        */
 428             i, j, m, k, teller = 0, slice; /* current slice number                           */
 429     real        nr_frames = 0;
 430     int*        slCount;      /* nr. of atoms in one slice                      */
 431     real        sdbangle = 0; /* sum of these angles                            */
 432     gmx_bool    use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
 433     rvec        direction, com;
 434     int         comsize, distsize;
 435     int *       comidx = nullptr, *distidx = nullptr;
 436     char*       grpname = nullptr;
 437     t_pbc       pbc;
 438     real        arcdist, tmpdist;
 439     gmx_rmpbc_t gpbc = nullptr;
 440
 441     /* PBC added for center-of-mass vector*/
 442     /* Initiate the pbc structure */
 443     std::memset(&pbc, 0, sizeof(pbc));
 444
 445     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 446     {
 447         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 448     }
 449
 450     nr_tails = index[1] - index[0];
 451     fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
 452     /* take first group as standard. Not rocksolid, but might catch error in index*/
 453
 454     if (permolecule)
 455     {
 456         nslices = nr_tails;
 457         bSliced = FALSE; /*force slices off */
 458         fprintf(stderr, "Calculating order parameters for each of %d molecules\n", nslices);
 459     }
 460
 461     if (radial)
 462     {
 463         use_unitvector = TRUE;
 464         fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
 465         get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
 466     }
 467     if (distcalc)
 468     {
 469         if (grpname != nullptr)
 470         {
 471             sfree(grpname);
 472         }
 473         fprintf(stderr, "Select an index group to use as distance reference\n");
 474         get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
 475         bSliced = FALSE; /*force slices off*/
 476     }
 477
 478     if (use_unitvector && bSliced)
 479     {
 480         fprintf(stderr,
 481                 "Warning:  slicing and specified unit vectors are not currently compatible\n");
 482     }
 483
 484     snew(slCount, nslices);
 485     snew(*slOrder, nslices);
 486     for (i = 0; i < nslices; i++)
 487     {
 488         snew((*slOrder)[i], ngrps);
 489     }
 490     if (distcalc)
 491     {
 492         snew(*distvals, nslices);
 493         for (i = 0; i < nslices; i++)
 494         {
 495             snew((*distvals)[i], ngrps);
 496         }
 497     }
 498     snew(*order, ngrps);
 499     snew(slFrameorder, nslices);
 500     snew(x1, natoms);
 501
 502     if (bSliced)
 503     {
 504         *slWidth = box[axis][axis] / static_cast<real>(nslices);
 505         fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", nslices, *slWidth);
 506     }
 507
 508
 509 #if 0
 510     nr_tails = index[1] - index[0];
 511     fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
 512     /* take first group as standard. Not rocksolid, but might catch error
 513        in index*/
 514 #endif
 515
 516     teller = 0;
 517
 518     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 519     /*********** Start processing trajectory ***********/
 520     do
 521     {
 522         if (bSliced)
 523         {
 524             *slWidth = box[axis][axis] / static_cast<real>(nslices);
 525         }
 526         teller++;
 527
 528         set_pbc(&pbc, ePBC, box);
 529         gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
 530
 531         /* Now loop over all groups. There are ngrps groups, the order parameter can
 532            be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
 533            atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
 534            course, in this case index[i+1] -index[i] has to be the same for all
 535            groups, namely the number of tails. i just runs over all atoms in a tail,
 536            so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
 537          */
 538
 539
 540         if (radial)
 541         {
 542             /*center-of-mass determination*/
 543             com[XX] = 0.0;
 544             com[YY] = 0.0;
 545             com[ZZ] = 0.0;
 546             for (j = 0; j < comsize; j++)
 547             {
 548                 rvec_inc(com, x1[comidx[j]]);
 549             }
 550             svmul(1.0 / comsize, com, com);
 551         }
 552         rvec displacementFromReference;
 553         if (distcalc)
 554         {
 555             rvec dref = { 0.0, 0.0, 0.0 };
 556             for (j = 0; j < distsize; j++)
 557             {
 558                 rvec_inc(dref, x1[distidx[j]]);
 559             }
 560             svmul(1.0 / distsize, dref, dref);
 561             if (radial)
 562             {
 563                 pbc_dx(&pbc, dref, com, displacementFromReference);
 564                 unitv(displacementFromReference, displacementFromReference);
 565             }
 566         }
 567
 568         for (i = 1; i < ngrps - 1; i++)
 569         {
 570             clear_rvec(frameorder);
 571
 572             size = index[i + 1] - index[i];
 573             if (size != nr_tails)
 574             {
 575                 gmx_fatal(FARGS,
 576                           "grp %d does not have same number of"
 577                           " elements as grp 1\n",
 578                           i);
 579             }
 580
 581             for (j = 0; j < size; j++)
 582             {
 583                 if (radial)
 584                 /*create unit vector*/
 585                 {
 586                     pbc_dx(&pbc, x1[a[index[i] + j]], com, direction);
 587                     unitv(direction, direction);
 588                     /*DEBUG*/
 589                     /*if (j==0)
 590                         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],
 591                             direction[0],direction[1],direction[2]);*/
 592                 }
 593
 594                 if (bUnsat)
 595                 {
 596                     rvec dist;
 597                     /* Using convention for unsaturated carbons */
 598                     /* first get Sz, the vector from Cn to Cn+1 */
 599                     rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], dist);
 600                     length = norm(dist);
 601                     check_length(length, a[index[i] + j], a[index[i + 1] + j]);
 602                     svmul(1.0 / length, dist, Sz);
 603
 604                     /* this is actually the cosine of the angle between the double bond
 605                        and axis, because Sz is normalized and the two other components of
 606                        the axis on the bilayer are zero */
 607                     if (use_unitvector)
 608                     {
 609                         sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
 610                     }
 611                     else
 612                     {
 613                         sdbangle += std::acos(Sz[axis]);
 614                     }
 615                 }
 616                 else
 617                 {
 618                     rvec dist;
 619                     /* get vector dist(Cn-1,Cn+1) for tail atoms */
 620                     rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i - 1] + j]], dist);
 621                     length = norm(dist); /* determine distance between two atoms */
 622                     check_length(length, a[index[i - 1] + j], a[index[i + 1] + j]);
 623
 624                     svmul(1.0 / length, dist, Sz);
 625                     /* Sz is now the molecular axis Sz, normalized and all that */
 626                 }
 627
 628                 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
 629                    we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
 630                 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], tmp1);
 631                 rvec_sub(x1[a[index[i - 1] + j]], x1[a[index[i] + j]], tmp2);
 632                 cprod(tmp1, tmp2, Sx);
 633                 svmul(1.0 / norm(Sx), Sx, Sx);
 634
 635                 /* now we can get Sy from the outer product of Sx and Sz   */
 636                 cprod(Sz, Sx, Sy);
 637                 svmul(1.0 / norm(Sy), Sy, Sy);
 638
 639                 /* the square of cosine of the angle between dist and the axis.
 640                    Using the innerproduct, but two of the three elements are zero
 641                    Determine the sum of the orderparameter of all atoms in group
 642                  */
 643                 if (use_unitvector)
 644                 {
 645                     cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
 646                     cossum[YY] = gmx::square(iprod(Sy, direction));
 647                     cossum[ZZ] = gmx::square(iprod(Sz, direction));
 648                 }
 649                 else
 650                 {
 651                     cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
 652                     cossum[YY] = gmx::square(Sy[axis]);
 653                     cossum[ZZ] = gmx::square(Sz[axis]);
 654                 }
 655
 656                 for (m = 0; m < DIM; m++)
 657                 {
 658                     frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
 659                 }
 660
 661                 if (bSliced)
 662                 {
 663                     /* get average coordinate in box length for slicing,
 664                        determine which slice atom is in, increase count for that
 665                        slice. slFrameorder and slOrder are reals, not
 666                        rvecs. Only the component [axis] of the order tensor is
 667                        kept, until I find it necessary to know the others too
 668                      */
 669
 670                     z1    = x1[a[index[i - 1] + j]][axis];
 671                     z2    = x1[a[index[i + 1] + j]][axis];
 672                     z_ave = 0.5 * (z1 + z2);
 673                     slice = static_cast<int>((static_cast<real>(nslices) * z_ave) / box[axis][axis]);
 674                     while (slice < 0)
 675                     {
 676                         slice += static_cast<real>(nslices);
 677                     }
 678                     slice = slice % nslices;
 679
 680                     slCount[slice]++; /* determine slice, increase count */
 681
 682                     slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
 683                 }
 684                 else if (permolecule)
 685                 {
 686                     /*  store per-molecule order parameter
 687                      *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 *
 688                      * iprod(cossum,direction) - 1); following is for Scd order: */
 689                     (*slOrder)[j][i] += -1
 690                                         * (1.0 / 3.0 * (3 * cossum[XX] - 1)
 691                                            + 1.0 / 3.0 * 0.5 * (3.0 * cossum[YY] - 1));
 692                 }
 693                 if (distcalc)
 694                 {
 695                     if (radial)
 696                     {
 697                         /* bin order parameter by arc distance from reference group*/
 698                         arcdist = gmx_angle(displacementFromReference, direction);
 699                         (*distvals)[j][i] += arcdist;
 700                     }
 701                     else if (i == 1)
 702                     {
 703                         /* Want minimum lateral distance to first group calculated */
 704                         tmpdist = trace(box); /* should be max value */
 705                         for (k = 0; k < distsize; k++)
 706                         {
 707                             rvec displacement;
 708                             pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i] + j]], displacement);
 709                             /* at the moment, just remove displacement[axis] */
 710                             displacement[axis] = 0;
 711                             tmpdist            = std::min(tmpdist, norm2(displacement));
 712                         }
 713                         // fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
 714                         (*distvals)[j][i] += std::sqrt(tmpdist);
 715                     }
 716                 }
 717             } /* end loop j, over all atoms in group */
 718
 719             for (m = 0; m < DIM; m++)
 720             {
 721                 (*order)[i][m] += (frameorder[m] / static_cast<real>(size));
 722             }
 723
 724             if (!permolecule)
 725             { /*Skip following if doing per-molecule*/
 726                 for (k = 0; k < nslices; k++)
 727                 {
 728                     if (slCount[k]) /* if no elements, nothing has to be added */
 729                     {
 730                         (*slOrder)[k][i] += slFrameorder[k] / static_cast<real>(slCount[k]);
 731                         slFrameorder[k] = 0;
 732                         slCount[k]      = 0;
 733                     }
 734                 }
 735             } /* end loop i, over all groups in indexfile */
 736         }
 737         nr_frames++;
 738
 739     } while (read_next_x(oenv, status, &t, x0, box));
 740     /*********** done with status file **********/
 741
 742     fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
 743     gmx_rmpbc_done(gpbc);
 744
 745     /* average over frames */
 746     for (i = 1; i < ngrps - 1; i++)
 747     {
 748         svmul(1.0 / nr_frames, (*order)[i], (*order)[i]);
 749         fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX], (*order)[i][YY],
 750                 (*order)[i][ZZ]);
 751         if (bSliced || permolecule)
 752         {
 753             for (k = 0; k < nslices; k++)
 754             {
 755                 (*slOrder)[k][i] /= nr_frames;
 756             }
 757         }
 758         if (distcalc)
 759         {
 760             for (k = 0; k < nslices; k++)
 761             {
 762                 (*distvals)[k][i] /= nr_frames;
 763             }
 764         }
 765     }
 766
 767     if (bUnsat)
 768     {
 769         fprintf(stderr, "Average angle between double bond and normal: %f\n",
 770                 180 * sdbangle / (nr_frames * static_cast<real>(size) * M_PI));
 771     }
 772
 773     sfree(x0); /* free memory used by coordinate arrays */
 774     sfree(x1);
 775     if (comidx != nullptr)
 776     {
 777         sfree(comidx);
 778     }
 779     if (distidx != nullptr)
 780     {
 781         sfree(distidx);
 782     }
 783     if (grpname != nullptr)
 784     {
 785         sfree(grpname);
 786     }
 787 }
 788
 789
 790 static void order_plot(rvec                    order[],
 791                        real*                   slOrder[],
 792                        const char*             afile,
 793                        const char*             bfile,
 794                        const char*             cfile,
 795                        int                     ngrps,
 796                        int                     nslices,
 797                        real                    slWidth,
 798                        gmx_bool                bSzonly,
 799                        gmx_bool                permolecule,
 800                        real**                  distvals,
 801                        const gmx_output_env_t* oenv)
 802 {
 803     FILE *ord, *slOrd; /* xvgr files with order parameters  */
 804     int   atom, slice; /* atom corresponding to order para.*/
 805     char  buf[256];    /* for xvgr title */
 806     real  S;           /* order parameter averaged over all atoms */
 807
 808     if (permolecule)
 809     {
 810         sprintf(buf, "Scd order parameters");
 811         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 812         sprintf(buf, "Orderparameters per atom per slice");
 813         slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
 814         for (atom = 1; atom < ngrps - 1; atom++)
 815         {
 816             fprintf(ord, "%12d   %12g\n", atom,
 817                     -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
 818         }
 819
 820         for (slice = 0; slice < nslices; slice++)
 821         {
 822             fprintf(slOrd, "%12d\t", slice);
 823             if (distvals)
 824             {
 825                 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
 826             }
 827             for (atom = 1; atom < ngrps - 1; atom++)
 828             {
 829                 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
 830             }
 831             fprintf(slOrd, "\n");
 832         }
 833     }
 834     else if (bSzonly)
 835     {
 836         sprintf(buf, "Orderparameters Sz per atom");
 837         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 838         fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
 839
 840         sprintf(buf, "Orderparameters per atom per slice");
 841         slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
 842
 843         for (atom = 1; atom < ngrps - 1; atom++)
 844         {
 845             fprintf(ord, "%12d       %12g\n", atom, order[atom][ZZ]);
 846         }
 847
 848         for (slice = 0; slice < nslices; slice++)
 849         {
 850             S = 0;
 851             for (atom = 1; atom < ngrps - 1; atom++)
 852             {
 853                 S += slOrder[slice][atom];
 854             }
 855             fprintf(slOrd, "%12g     %12g\n", static_cast<real>(slice) * slWidth,
 856                     S / static_cast<real>(atom));
 857         }
 858     }
 859     else
 860     {
 861         sprintf(buf, "Order tensor diagonal elements");
 862         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 863         sprintf(buf, "Deuterium order parameters");
 864         slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
 865
 866         for (atom = 1; atom < ngrps - 1; atom++)
 867         {
 868             fprintf(ord, "%12d   %12g   %12g   %12g\n", atom, order[atom][XX], order[atom][YY],
 869                     order[atom][ZZ]);
 870             fprintf(slOrd, "%12d   %12g\n", atom,
 871                     -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
 872         }
 873     }
 874     xvgrclose(ord);
 875     xvgrclose(slOrd);
 876 }
 877
 878 static void write_bfactors(t_filenm*         fnm,
 879                            int               nfile,
 880                            const int*        index,
 881                            const int*        a,
 882                            int               nslices,
 883                            int               ngrps,
 884                            real**            order,
 885                            const t_topology* top,
 886                            real**            distvals,
 887                            gmx_output_env_t* oenv)
 888 {
 889     /*function to write order parameters as B factors in PDB file using
 890           first frame of trajectory*/
 891     t_trxstatus* status;
 892     t_trxframe   fr, frout;
 893     t_atoms      useatoms;
 894     int          i, j, ctr, nout;
 895
 896     ngrps -= 2; /*we don't have an order parameter for the first or
 897                       last atom in each chain*/
 898     nout = nslices * ngrps;
 899     read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
 900
 901     close_trx(status);
 902     frout        = fr;
 903     frout.natoms = nout;
 904     frout.bF     = FALSE;
 905     frout.bV     = FALSE;
 906     frout.x      = nullptr;
 907     snew(frout.x, nout);
 908
 909     init_t_atoms(&useatoms, nout, TRUE);
 910     useatoms.nr = nout;
 911
 912     /*initialize PDBinfo*/
 913     for (i = 0; i < useatoms.nr; ++i)
 914     {
 915         useatoms.pdbinfo[i].type         = 0;
 916         useatoms.pdbinfo[i].occup        = 0.0;
 917         useatoms.pdbinfo[i].bfac         = 0.0;
 918         useatoms.pdbinfo[i].bAnisotropic = FALSE;
 919     }
 920
 921     for (j = 0, ctr = 0; j < nslices; j++)
 922     {
 923         for (i = 0; i < ngrps; i++, ctr++)
 924         {
 925             /*iterate along each chain*/
 926             useatoms.pdbinfo[ctr].bfac = order[j][i + 1];
 927             if (distvals)
 928             {
 929                 useatoms.pdbinfo[ctr].occup = distvals[j][i + 1];
 930             }
 931             copy_rvec(fr.x[a[index[i + 1] + j]], frout.x[ctr]);
 932             useatoms.atomname[ctr] = top->atoms.atomname[a[index[i + 1] + j]];
 933             useatoms.atom[ctr]     = top->atoms.atom[a[index[i + 1] + j]];
 934             useatoms.nres          = std::max(useatoms.nres, useatoms.atom[ctr].resind + 1);
 935             useatoms.resinfo[useatoms.atom[ctr].resind] =
 936                     top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
 937         }
 938     }
 939
 940     write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, nullptr,
 941                    frout.ePBC, frout.box);
 942
 943     sfree(frout.x);
 944     done_atom(&useatoms);
 945 }
 946
 947 int gmx_order(int argc, char* argv[])
 948 {
 949     const char* desc[] = {
 950         "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
 951         "vector i-1, i+1 is used together with an axis. ",
 952         "The index file should contain only the groups to be used for calculations,",
 953         "with each group of equivalent carbons along the relevant acyl chain in its own",
 954         "group. There should not be any generic groups (like System, Protein) in the index",
 955         "file to avoid confusing the program (this is not relevant to tetrahedral order",
 956         "parameters however, which only work for water anyway).[PAR]",
 957         "[THISMODULE] can also give all",
 958         "diagonal elements of the order tensor and even calculate the deuterium",
 959         "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
 960         "order tensor component (specified by the [TT]-d[tt] option) is given and the",
 961         "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
 962         "selected, all diagonal elements and the deuterium order parameter is",
 963         "given.[PAR]",
 964         "The tetrahedrality order parameters can be determined",
 965         "around an atom. Both angle an distance order parameters are calculated. See",
 966         "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
 967         "for more details."
 968     };
 969
 970     static int         nslices       = 1;     /* nr of slices defined       */
 971     static gmx_bool    bSzonly       = FALSE; /* True if only Sz is wanted  */
 972     static gmx_bool    bUnsat        = FALSE; /* True if carbons are unsat. */
 973     static const char* normal_axis[] = { nullptr, "z", "x", "y", nullptr };
 974     static gmx_bool    permolecule   = FALSE; /*compute on a per-molecule basis */
 975     static gmx_bool    radial        = FALSE; /*compute a radial membrane normal */
 976     static gmx_bool    distcalc      = FALSE; /*calculate distance from a reference group */
 977     t_pargs            pa[]          = {
 978         { "-d", FALSE, etENUM, { normal_axis }, "Direction of the normal on the membrane" },
 979         { "-sl",
 980           FALSE,
 981           etINT,
 982           { &nslices },
 983           "Calculate order parameter as function of box length, dividing the box"
 984           " into this number of slices." },
 985         { "-szonly",
 986           FALSE,
 987           etBOOL,
 988           { &bSzonly },
 989           "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
 990         { "-unsat",
 991           FALSE,
 992           etBOOL,
 993           { &bUnsat },
 994           "Calculate order parameters for unsaturated carbons. Note that this can"
 995           "not be mixed with normal order parameters." },
 996         { "-permolecule",
 997           FALSE,
 998           etBOOL,
 999           { &permolecule },
1000           "Compute per-molecule Scd order parameters" },
1001         { "-radial", FALSE, etBOOL, { &radial }, "Compute a radial membrane normal" },
1002         { "-calcdist", FALSE, etBOOL, { &distcalc }, "Compute distance from a reference" },
1003     };
1004
1005     rvec*  order;         /* order par. for each atom   */
1006     real** slOrder;       /* same, per slice            */
1007     real   slWidth = 0.0; /* width of a slice           */
1008     char** grpname;       /* groupnames                 */
1009     int    ngrps,         /* nr. of groups              */
1010             i, axis = 0;  /* normal axis                */
1011     t_topology* top;      /* topology         */
1012     int         ePBC;
1013     int *       index, /* indices for a              */
1014             *a;        /* atom numbers in each group */
1015     t_blocka* block;   /* data from index file       */
1016     t_filenm  fnm[] = {
1017         /* files for g_order    */
1018         { efTRX, "-f", nullptr, ffREAD },     /* trajectory file              */
1019         { efNDX, "-n", nullptr, ffREAD },     /* index file           */
1020         { efNDX, "-nr", nullptr, ffOPTRD },   /* index for radial axis calculation */
1021         { efTPR, nullptr, nullptr, ffREAD },  /* topology file                */
1022         { efXVG, "-o", "order", ffWRITE },    /* xvgr output file     */
1023         { efXVG, "-od", "deuter", ffWRITE },  /* xvgr output file           */
1024         { efPDB, "-ob", nullptr, ffOPTWR },   /* write Scd as B factors to PDB if permolecule   */
1025         { efXVG, "-os", "sliced", ffWRITE },  /* xvgr output file           */
1026         { efXVG, "-Sg", "sg-ang", ffOPTWR },  /* xvgr output file           */
1027         { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file           */
1028         { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR },  /* xvgr output file           */
1029         { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file           */
1030     };
1031     gmx_bool bSliced = FALSE; /* True if box is sliced      */
1032 #define NFILE asize(fnm)
1033     real**            distvals = nullptr;
1034     const char *      sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
1035     gmx_output_env_t* oenv;
1036
1037     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
1038                            asize(desc), desc, 0, nullptr, &oenv))
1039     {
1040         return 0;
1041     }
1042     if (nslices < 1)
1043     {
1044         gmx_fatal(FARGS, "Can not have nslices < 1");
1045     }
1046     sgfnm  = opt2fn_null("-Sg", NFILE, fnm);
1047     skfnm  = opt2fn_null("-Sk", NFILE, fnm);
1048     ndxfnm = opt2fn_null("-n", NFILE, fnm);
1049     tpsfnm = ftp2fn(efTPR, NFILE, fnm);
1050     trxfnm = ftp2fn(efTRX, NFILE, fnm);
1051
1052     /* Calculate axis */
1053     GMX_RELEASE_ASSERT(normal_axis[0] != nullptr, "Options inconsistency; normal_axis[0] is NULL");
1054     if (std::strcmp(normal_axis[0], "x") == 0)
1055     {
1056         axis = XX;
1057     }
1058     else if (std::strcmp(normal_axis[0], "y") == 0)
1059     {
1060         axis = YY;
1061     }
1062     else if (std::strcmp(normal_axis[0], "z") == 0)
1063     {
1064         axis = ZZ;
1065     }
1066     else
1067     {
1068         gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1069     }
1070
1071     switch (axis)
1072     {
1073         case 0: fprintf(stderr, "Taking x axis as normal to the membrane\n"); break;
1074         case 1: fprintf(stderr, "Taking y axis as normal to the membrane\n"); break;
1075         case 2: fprintf(stderr, "Taking z axis as normal to the membrane\n"); break;
1076     }
1077
1078     /* tetraheder order parameter */
1079     if (skfnm || sgfnm)
1080     {
1081         /* If either of theoptions is set we compute both */
1082         sgfnm = opt2fn("-Sg", NFILE, fnm);
1083         skfnm = opt2fn("-Sk", NFILE, fnm);
1084         calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1085                               opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm), oenv);
1086         /* view xvgr files */
1087         do_view(oenv, opt2fn("-Sg", NFILE, fnm), nullptr);
1088         do_view(oenv, opt2fn("-Sk", NFILE, fnm), nullptr);
1089         if (nslices > 1)
1090         {
1091             do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), nullptr);
1092             do_view(oenv, opt2fn("-Sksl", NFILE, fnm), nullptr);
1093         }
1094     }
1095     else
1096     {
1097         /* tail order parameter */
1098
1099         if (nslices > 1)
1100         {
1101             bSliced = TRUE;
1102             fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1103         }
1104
1105         if (bSzonly)
1106         {
1107             fprintf(stderr, "Only calculating Sz\n");
1108         }
1109         if (bUnsat)
1110         {
1111             fprintf(stderr, "Taking carbons as unsaturated!\n");
1112         }
1113
1114         top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
1115
1116         block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1117         index = block->index; /* get indices from t_block block */
1118         a     = block->a;     /* see block.h                    */
1119         ngrps = block->nr;
1120
1121         if (permolecule)
1122         {
1123             nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1124             fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1125         }
1126
1127         if (radial)
1128         {
1129             fprintf(stderr, "Calculating radial distances\n");
1130             if (!permolecule)
1131             {
1132                 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1133             }
1134         }
1135
1136         /* show atomtypes, to check if index file is correct */
1137         print_types(index, a, ngrps, grpname, top);
1138
1139         calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order, &slOrder, &slWidth, nslices,
1140                    bSliced, bUnsat, top, ePBC, ngrps, axis, permolecule, radial, distcalc,
1141                    opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
1142
1143         if (radial)
1144         {
1145             ngrps--; /*don't print the last group--was used for
1146                                center-of-mass determination*/
1147         }
1148         order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
1149                    opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule,
1150                    distvals, oenv);
1151
1152         if (opt2bSet("-ob", NFILE, fnm))
1153         {
1154             if (!permolecule)
1155             {
1156                 fprintf(stderr,
1157                         "Won't write B-factors with averaged order parameters; use -permolecule\n");
1158             }
1159             else
1160             {
1161                 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1162             }
1163         }
1164
1165
1166         do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);  /* view xvgr file */
1167         do_view(oenv, opt2fn("-os", NFILE, fnm), nullptr); /* view xvgr file */
1168         do_view(oenv, opt2fn("-od", NFILE, fnm), nullptr); /* view xvgr file */
1169     }
1170
1171     if (distvals != nullptr)
1172     {
1173         for (i = 0; i < nslices; ++i)
1174         {
1175             sfree(distvals[i]);
1176         }
1177         sfree(distvals);
1178     }
1179
1180     return 0;
1181 }