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