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, 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, 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 = NULL, *fpsk = NULL;
 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 = NULL;
 275
 276
 277     read_tps_conf(fnTPS, &top, &ePBC, &xtop, NULL, 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(NULL, ng, index[0], NULL, 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_trj(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(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 void calc_order(const char *fn, 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  = NULL, *distidx = NULL;
 409     char        *grpname = NULL;
 410     t_pbc        pbc;
 411     real         arcdist, tmpdist;
 412     gmx_rmpbc_t  gpbc = NULL;
 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 != NULL)
 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                     if (z_ave < 0)
 641                     {
 642                         z_ave += box[axis][axis];
 643                     }
 644                     if (z_ave > box[axis][axis])
 645                     {
 646                         z_ave -= box[axis][axis];
 647                     }
 648
 649                     slice  = static_cast<int>((0.5 + (z_ave / (*slWidth))) - 1);
 650                     slCount[slice]++;     /* determine slice, increase count */
 651
 652                     slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
 653                 }
 654                 else if (permolecule)
 655                 {
 656                     /*  store per-molecule order parameter
 657                      *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
 658                      *  following is for Scd order: */
 659                     (*slOrder)[j][i] += -1* (1.0/3.0 * (3 * cossum[XX] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum[YY] - 1));
 660                 }
 661                 if (distcalc)
 662                 {
 663                     if (radial)
 664                     {
 665                         /* bin order parameter by arc distance from reference group*/
 666                         arcdist            = gmx_angle(dvec, direction);
 667                         (*distvals)[j][i] += arcdist;
 668                     }
 669                     else if (i == 1)
 670                     {
 671                         /* Want minimum lateral distance to first group calculated */
 672                         tmpdist = trace(box);  /* should be max value */
 673                         for (k = 0; k < distsize; k++)
 674                         {
 675                             pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
 676                             /* at the moment, just remove dvec[axis] */
 677                             dvec[axis] = 0;
 678                             tmpdist    = std::min(tmpdist, norm2(dvec));
 679                         }
 680                         //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
 681                         (*distvals)[j][i] += std::sqrt(tmpdist);
 682                     }
 683                 }
 684             } /* end loop j, over all atoms in group */
 685
 686             for (m = 0; m < DIM; m++)
 687             {
 688                 (*order)[i][m] += (frameorder[m]/size);
 689             }
 690
 691             if (!permolecule)
 692             {   /*Skip following if doing per-molecule*/
 693                 for (k = 0; k < nslices; k++)
 694                 {
 695                     if (slCount[k]) /* if no elements, nothing has to be added */
 696                     {
 697                         (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
 698                         slFrameorder[k]   = 0; slCount[k] = 0;
 699                     }
 700                 }
 701             } /* end loop i, over all groups in indexfile */
 702         }
 703         nr_frames++;
 704
 705     }
 706     while (read_next_x(oenv, status, &t, x0, box));
 707     /*********** done with status file **********/
 708
 709     fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
 710     gmx_rmpbc_done(gpbc);
 711
 712     /* average over frames */
 713     for (i = 1; i < ngrps - 1; i++)
 714     {
 715         svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
 716         fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
 717                 (*order)[i][YY], (*order)[i][ZZ]);
 718         if (bSliced || permolecule)
 719         {
 720             for (k = 0; k < nslices; k++)
 721             {
 722                 (*slOrder)[k][i] /= nr_frames;
 723             }
 724         }
 725         if (distcalc)
 726         {
 727             for (k = 0; k < nslices; k++)
 728             {
 729                 (*distvals)[k][i] /= nr_frames;
 730             }
 731         }
 732     }
 733
 734     if (bUnsat)
 735     {
 736         fprintf(stderr, "Average angle between double bond and normal: %f\n",
 737                 180*sdbangle/(nr_frames * size*M_PI));
 738     }
 739
 740     sfree(x0); /* free memory used by coordinate arrays */
 741     sfree(x1);
 742     if (comidx != NULL)
 743     {
 744         sfree(comidx);
 745     }
 746     if (distidx != NULL)
 747     {
 748         sfree(distidx);
 749     }
 750     if (grpname != NULL)
 751     {
 752         sfree(grpname);
 753     }
 754 }
 755
 756
 757 void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
 758                 const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
 759                 gmx_bool permolecule, real **distvals, const gmx_output_env_t *oenv)
 760 {
 761     FILE       *ord, *slOrd;      /* xvgr files with order parameters  */
 762     int         atom, slice;      /* atom corresponding to order para.*/
 763     char        buf[256];         /* for xvgr title */
 764     real        S;                /* order parameter averaged over all atoms */
 765
 766     if (permolecule)
 767     {
 768         sprintf(buf, "Scd order parameters");
 769         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 770         sprintf(buf, "Orderparameters per atom per slice");
 771         slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
 772         for (atom = 1; atom < ngrps - 1; atom++)
 773         {
 774             fprintf(ord, "%12d   %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
 775                                                         1.0/3.0 * order[atom][YY]));
 776         }
 777
 778         for (slice = 0; slice < nslices; slice++)
 779         {
 780             fprintf(slOrd, "%12d\t", slice);
 781             if (distvals)
 782             {
 783                 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
 784             }
 785             for (atom = 1; atom < ngrps - 1; atom++)
 786             {
 787                 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
 788             }
 789             fprintf(slOrd, "\n");
 790         }
 791
 792     }
 793     else if (bSzonly)
 794     {
 795         sprintf(buf, "Orderparameters Sz per atom");
 796         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 797         fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
 798
 799         sprintf(buf, "Orderparameters per atom per slice");
 800         slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
 801
 802         for (atom = 1; atom < ngrps - 1; atom++)
 803         {
 804             fprintf(ord, "%12d       %12g\n", atom, order[atom][ZZ]);
 805         }
 806
 807         for (slice = 0; slice < nslices; slice++)
 808         {
 809             S = 0;
 810             for (atom = 1; atom < ngrps - 1; atom++)
 811             {
 812                 S += slOrder[slice][atom];
 813             }
 814             fprintf(slOrd, "%12g     %12g\n", slice*slWidth, S/atom);
 815         }
 816
 817     }
 818     else
 819     {
 820         sprintf(buf, "Order tensor diagonal elements");
 821         ord = xvgropen(afile, buf, "Atom", "S", oenv);
 822         sprintf(buf, "Deuterium order parameters");
 823         slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
 824
 825         for (atom = 1; atom < ngrps - 1; atom++)
 826         {
 827             fprintf(ord, "%12d   %12g   %12g   %12g\n", atom, order[atom][XX],
 828                     order[atom][YY], order[atom][ZZ]);
 829             fprintf(slOrd, "%12d   %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
 830                                                           1.0/3.0 * order[atom][YY]));
 831         }
 832
 833     }
 834     xvgrclose(ord);
 835     xvgrclose(slOrd);
 836 }
 837
 838 void write_bfactors(t_filenm  *fnm, int nfile, int *index, int *a, int nslices, int ngrps, real **order, const t_topology *top, real **distvals, gmx_output_env_t *oenv)
 839 {
 840     /*function to write order parameters as B factors in PDB file using
 841           first frame of trajectory*/
 842     t_trxstatus *status;
 843     t_trxframe   fr, frout;
 844     t_atoms      useatoms;
 845     int          i, j, ctr, nout;
 846
 847     ngrps -= 2;  /*we don't have an order parameter for the first or
 848                        last atom in each chain*/
 849     nout   = nslices*ngrps;
 850     read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
 851
 852     close_trj(status);
 853     frout        = fr;
 854     frout.natoms = nout;
 855     frout.bF     = FALSE;
 856     frout.bV     = FALSE;
 857     frout.x      = 0;
 858     snew(frout.x, nout);
 859
 860     init_t_atoms(&useatoms, nout, TRUE);
 861     useatoms.nr = nout;
 862
 863     /*initialize PDBinfo*/
 864     for (i = 0; i < useatoms.nr; ++i)
 865     {
 866         useatoms.pdbinfo[i].type         = 0;
 867         useatoms.pdbinfo[i].occup        = 0.0;
 868         useatoms.pdbinfo[i].bfac         = 0.0;
 869         useatoms.pdbinfo[i].bAnisotropic = FALSE;
 870     }
 871
 872     for (j = 0, ctr = 0; j < nslices; j++)
 873     {
 874         for (i = 0; i < ngrps; i++, ctr++)
 875         {
 876             /*iterate along each chain*/
 877             useatoms.pdbinfo[ctr].bfac = order[j][i+1];
 878             if (distvals)
 879             {
 880                 useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
 881             }
 882             copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
 883             useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
 884             useatoms.atom[ctr]     = top->atoms.atom[a[index[i+1]+j]];
 885             useatoms.nres          = std::max(useatoms.nres, useatoms.atom[ctr].resind+1);
 886             useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
 887         }
 888     }
 889
 890     write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);
 891
 892     sfree(frout.x);
 893     done_atom(&useatoms);
 894 }
 895
 896 int gmx_order(int argc, char *argv[])
 897 {
 898     const char        *desc[] = {
 899         "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
 900         "vector i-1, i+1 is used together with an axis. ",
 901         "The index file should contain only the groups to be used for calculations,",
 902         "with each group of equivalent carbons along the relevant acyl chain in its own",
 903         "group. There should not be any generic groups (like System, Protein) in the index",
 904         "file to avoid confusing the program (this is not relevant to tetrahedral order",
 905         "parameters however, which only work for water anyway).[PAR]",
 906         "[THISMODULE] can also give all",
 907         "diagonal elements of the order tensor and even calculate the deuterium",
 908         "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
 909         "order tensor component (specified by the [TT]-d[tt] option) is given and the",
 910         "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
 911         "selected, all diagonal elements and the deuterium order parameter is",
 912         "given.[PAR]"
 913         "The tetrahedrality order parameters can be determined",
 914         "around an atom. Both angle an distance order parameters are calculated. See",
 915         "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
 916         "for more details."
 917     };
 918
 919     static int         nslices       = 1;     /* nr of slices defined       */
 920     static gmx_bool    bSzonly       = FALSE; /* True if only Sz is wanted  */
 921     static gmx_bool    bUnsat        = FALSE; /* True if carbons are unsat. */
 922     static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
 923     static gmx_bool    permolecule   = FALSE; /*compute on a per-molecule basis */
 924     static gmx_bool    radial        = FALSE; /*compute a radial membrane normal */
 925     static gmx_bool    distcalc      = FALSE; /*calculate distance from a reference group */
 926     t_pargs            pa[]          = {
 927         { "-d",      FALSE, etENUM, {normal_axis},
 928           "Direction of the normal on the membrane" },
 929         { "-sl",     FALSE, etINT, {&nslices},
 930           "Calculate order parameter as function of box length, dividing the box"
 931           " into this number of slices." },
 932         { "-szonly", FALSE, etBOOL, {&bSzonly},
 933           "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
 934         { "-unsat",  FALSE, etBOOL, {&bUnsat},
 935           "Calculate order parameters for unsaturated carbons. Note that this can"
 936           "not be mixed with normal order parameters." },
 937         { "-permolecule", FALSE, etBOOL, {&permolecule},
 938           "Compute per-molecule Scd order parameters" },
 939         { "-radial", FALSE, etBOOL, {&radial},
 940           "Compute a radial membrane normal" },
 941         { "-calcdist", FALSE, etBOOL, {&distcalc},
 942           "Compute distance from a reference" },
 943     };
 944
 945     rvec              *order;                         /* order par. for each atom   */
 946     real             **slOrder;                       /* same, per slice            */
 947     real               slWidth = 0.0;                 /* width of a slice           */
 948     char             **grpname;                       /* groupnames                 */
 949     int                ngrps,                         /* nr. of groups              */
 950                        i,
 951                        axis = 0;                      /* normal axis                */
 952     t_topology       *top;                            /* topology         */
 953     int               ePBC;
 954     int              *index,                          /* indices for a              */
 955     *a;                                               /* atom numbers in each group */
 956     t_blocka         *block;                          /* data from index file       */
 957     t_filenm          fnm[] = {                       /* files for g_order    */
 958         { efTRX, "-f", NULL,  ffREAD },               /* trajectory file              */
 959         { efNDX, "-n", NULL,  ffREAD },               /* index file           */
 960         { efNDX, "-nr", NULL,  ffREAD },              /* index for radial axis calculation        */
 961         { efTPR, NULL, NULL,  ffREAD },               /* topology file                */
 962         { efXVG, "-o", "order", ffWRITE },            /* xvgr output file     */
 963         { efXVG, "-od", "deuter", ffWRITE },          /* xvgr output file           */
 964         { efPDB, "-ob", NULL, ffWRITE },              /* write Scd as B factors to PDB if permolecule           */
 965         { efXVG, "-os", "sliced", ffWRITE },          /* xvgr output file           */
 966         { efXVG, "-Sg", "sg-ang", ffOPTWR },          /* xvgr output file           */
 967         { efXVG, "-Sk", "sk-dist", ffOPTWR },         /* xvgr output file           */
 968         { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR },  /* xvgr output file           */
 969         { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file           */
 970     };
 971     gmx_bool          bSliced = FALSE;                /* True if box is sliced      */
 972 #define NFILE asize(fnm)
 973     real            **distvals = NULL;
 974     const char       *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
 975     gmx_output_env_t *oenv;
 976
 977     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 978                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 979     {
 980         return 0;
 981     }
 982     if (nslices < 1)
 983     {
 984         gmx_fatal(FARGS, "Can not have nslices < 1");
 985     }
 986     sgfnm  = opt2fn_null("-Sg", NFILE, fnm);
 987     skfnm  = opt2fn_null("-Sk", NFILE, fnm);
 988     ndxfnm = opt2fn_null("-n", NFILE, fnm);
 989     tpsfnm = ftp2fn(efTPR, NFILE, fnm);
 990     trxfnm = ftp2fn(efTRX, NFILE, fnm);
 991
 992     /* Calculate axis */
 993     GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Options inconsistency; normal_axis[0] is NULL");
 994     if (std::strcmp(normal_axis[0], "x") == 0)
 995     {
 996         axis = XX;
 997     }
 998     else if (std::strcmp(normal_axis[0], "y") == 0)
 999     {
1000         axis = YY;
1001     }
1002     else if (std::strcmp(normal_axis[0], "z") == 0)
1003     {
1004         axis = ZZ;
1005     }
1006     else
1007     {
1008         gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1009     }
1010
1011     switch (axis)
1012     {
1013         case 0:
1014             fprintf(stderr, "Taking x axis as normal to the membrane\n");
1015             break;
1016         case 1:
1017             fprintf(stderr, "Taking y axis as normal to the membrane\n");
1018             break;
1019         case 2:
1020             fprintf(stderr, "Taking z axis as normal to the membrane\n");
1021             break;
1022     }
1023
1024     /* tetraheder order parameter */
1025     if (skfnm || sgfnm)
1026     {
1027         /* If either of theoptions is set we compute both */
1028         sgfnm = opt2fn("-Sg", NFILE, fnm);
1029         skfnm = opt2fn("-Sk", NFILE, fnm);
1030         calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1031                               opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm),
1032                               oenv);
1033         /* view xvgr files */
1034         do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
1035         do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
1036         if (nslices > 1)
1037         {
1038             do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
1039             do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
1040         }
1041     }
1042     else
1043     {
1044         /* tail order parameter */
1045
1046         if (nslices > 1)
1047         {
1048             bSliced = TRUE;
1049             fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1050         }
1051
1052         if (bSzonly)
1053         {
1054             fprintf(stderr, "Only calculating Sz\n");
1055         }
1056         if (bUnsat)
1057         {
1058             fprintf(stderr, "Taking carbons as unsaturated!\n");
1059         }
1060
1061         top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
1062
1063         block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1064         index = block->index;                   /* get indices from t_block block */
1065         a     = block->a;                       /* see block.h                    */
1066         ngrps = block->nr;
1067
1068         if (permolecule)
1069         {
1070             nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1071             fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1072         }
1073
1074         if (radial)
1075         {
1076             fprintf(stderr, "Calculating radial distances\n");
1077             if (!permolecule)
1078             {
1079                 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1080             }
1081         }
1082
1083         /* show atomtypes, to check if index file is correct */
1084         print_types(index, a, ngrps, grpname, top);
1085
1086         calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order,
1087                    &slOrder, &slWidth, nslices, bSliced, bUnsat,
1088                    top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
1089
1090         if (radial)
1091         {
1092             ngrps--; /*don't print the last group--was used for
1093                                center-of-mass determination*/
1094
1095         }
1096         order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
1097                    opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);
1098
1099         if (opt2bSet("-ob", NFILE, fnm))
1100         {
1101             if (!permolecule)
1102             {
1103                 fprintf(stderr,
1104                         "Won't write B-factors with averaged order parameters; use -permolecule\n");
1105             }
1106             else
1107             {
1108                 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1109             }
1110         }
1111
1112
1113         do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);  /* view xvgr file */
1114         do_view(oenv, opt2fn("-os", NFILE, fnm), NULL); /* view xvgr file */
1115         do_view(oenv, opt2fn("-od", NFILE, fnm), NULL); /* view xvgr file */
1116     }
1117
1118     if (distvals != NULL)
1119     {
1120         for (i = 0; i < nslices; ++i)
1121         {
1122             sfree(distvals[i]);
1123         }
1124         sfree(distvals);
1125     }
1126
1127     return 0;
1128 }