src/gromacs/gmxana/gmx_order.c

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