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